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In  this  dissertation  we  discuss  the  numerical  solution  of  systems  of  hyperbolic  partial 
differential  equations  with  lower  order  terms  and  step  function  initial  data.  These  equations 
arise  in  modeling  the  propagation  of  a  signal  with  loss,  such  as  a  signal  in  a  resistive 
co-axial  cable,  or  the  flow  of  neutrons  in  a  reactor.  Majda  and  Osher  have  shown  that 
dissipative  finite  difference  approximations  to  such  problems  display  a  numerical  artifact 
which  is  not  encountered  for  scalar  equations.  Namely,  noise  from  an  initial  discontinuity 
propagates  into  a  large  region  behind  the  discontinuity.  Their  results  do  not  apply  in 
the  vicinity  of  a  discontinuity,  and  our  goal  is  to  discover  the  detailed  behavior  in  this 
region.  This  information  will  be  of  use  in  constructing  algorithms  that  attempt  to  accurately 
approximate  solutions  with  discontinuities  or  shocks.  . 

We  analyze  the  behavior  of  finite  difference  approximations  to  a  particular  model 
problem  with  step  function  initial  data.  We  use  Fourier  transforms  to  express  the  solution 
of  the  difference  approximation  as  a  Fourier  integral.  The  behavior  of  this  integral  is  then 
examined  by  using  the  theory  of  uniform  asymptotic  estimation  of  integrals.  We  show  that 
the  solution  of  the  difference  approximation  is  modeled  by  a  particular  class  of  integrals, 
which  we  call  generalized  Bessel  functions.  With  these  results,  we  are  able  to  show  that 
there  is  a  “smearing”  of  the  discontinuity  which  is  the  same  as  that  which  one  obtains 
when  approximating  the  scalar  wave  equation.  Thus,  the  behavior  near  the  discontinuity  is 
qualitatively  similiar  to  the  near-front  behavior  of  the  scalar  wave  equation.  These  results 
show  that  the  artifact  discovered  by  Majda  and  Osher  is  overwhelmed  near  the  discontinuity 
by  the  numerical  dispersion  and/or  dissipation  introduced  by  the  difference  approximation. 
Graphs  of  the  generalized  Bessel  functions  are  provided.  Finally,  we  discuss  the  relevance 
of  our  results  to  automatic  mesh  refinement  and  give  an  example. 
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1.  Introduction 


In  this  dissertation  we  consider  the  numerical  solution  of  the  telegrapher’s  equation 


du  dv 

dt  dx 

dv  du 

dt  dx 


(11) 


with  step  function  initial  data.  This  equation  models  wave  propagation  in  a  system  of 
equations  with  coupling  caused  by  the  presence  of  undifferentiated  or  lower  order  terms. 
(If  the  lower  order  term  were  not  present,  it  would  be  possible  to  diagonalize  the  system 
with  a  linear  change  of  variables.)  This  equation  is  of  interest  because  the  solution  to 
this  equation  is  smooth  away  from  discontinuities  in  the  solution  (Birkhoff  and  Lynch 
[1966])  but  numerical  approximations  are  not  (Majda  and  Osher  [1977]).  Also,  as  we  show 
below,  the  equation  itself  can  be  considered  a  model  for  transport  phenomena,  so  that  a 
better  understanding  of  the  behavior  of  approximate  solutions  to  it  will  help  in  interpreting 
numerical  results  and  in  constructing  numerical  methods. 

In  this  section  we  will  first  discuss  some  examples  from  engineering  and  physics  of 
equations  with  coupling  through  lower  order  terms.  These  examples  will  help  develop  an 
intuition  about  what  the  exact  solutions  should  look  like.  We  will  then  give  a  brief  historical 
overview  of  material  on  both  the  telegrapher’s  equation  and  on  techniques  for  analyzing 
the  behavior  of  approximations  to  hyperbolic  partial  differential  equations  (to  which  class 
the  telegrapher’s  equation  belongs).  Finally,  we  will  discuss  in  more  detail  the  methods 
which  we  shall  use  to  find  the  behavior  of  approximations  to  (1.1)  near  the  discontinuity. 

The  following  is  an  overview  of  this  dissertation.  The  theoretical  results  are  presented 
in  sections  2,  3,  and  4.  In  section  2,  we  find  an  asymptotic  representation  to  the  solution  of 
(1.1)  near  the  step  discontinuity.  The  purpose  of  this  section  is  to  both  show  the  behavior 
of  (1.1)  to  which  the  approximate  solutions  will  be  compared  and  to  illustrate  the  methods 
which  will  be  used  in  later  sections  to  analyze  difference  approximations  to  (1.1).  In  section 
3,  we  analyze  a  semidiscrete  second  order  centered  difference  approximation  to  (1.1).  We 
first  show  that  the  width  of  the  front  in  the  approximation  is  0(h2/3).  We  then  derive 
a  form  of  integral  (a  generalized  Bessel  function)  which  represents  the  solution  near  the 
front.  Finally,  we  compare  the  behavior  of  the  approximation  to  (1.1)  with  the  same  type 
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of  approximation  to  the  scalar  wave  equation.  In  section  4,  we  analyze  a  model  equation 
which  represents  a  general  difference  approximation  to  (1.1).  This  section  generalizes 
the  results  from  section  3.  Section  5  contains  numerical  results.  We  consider  here  four 
different  difference  approximations  to  (1.1).  Graphs  of  the  computed  solutions  and  the 
observed  widths  of  the  front  for  each  method  are  presented.  These  are  followed  by  graphs 
of  the  generalized  Bessel  functions  which  we  show  in  section  4  to  represent  the  solution 
to  a  difference  approximation  to  (1.1)  near  the  front.  Finally,  an  application  of  the  results 
is  made  to  a  computation  using  mesh  refinement. 


Examples  of  Equations  with  Lower  Order  Terms 

We  begin  with  a  brief  discussion  of  a  few  examples  of  equations  with  lower  order 
terms  to  show  where  the  problem  we  address  arises.  The  basic  physical  feature  of  all  of 
these  equations  is  a  finite  speed  of  propagation  combined  with  decay  of  signal  through, 
for  example,  absorption. 

Perhaps  the  most  familiar  example  is  transmission  along  a  resistive  coaxial  cable.  This 
equation  is  usually  written  as 

dzip  fit  d2ip  Azfio  dip 
a?  —  ^2  ~M2  +  dt 


(see  Morse  and  Feshbach  [1953],  pp.  218-9).  Here,  fi,  c,  and  tr  are  physical  constants  of 
the  cable,  and  c  is  the  speed  of  light,  ip  is  a  potential  for  the  fields,  /  is  time,  and  z  is 
length  along  the  cable.  This  equation  may  be  rewritten  as 


dip 

~di 

fit  d<p 
c2  dt 


dtp 

dz 

dip 

dz 


AltfltJ 


<p 


where  <p  is  introduced  only  to  convert  the  single  second  order  equation  into  two  first  order 
equations.  Another  example  is  transmission  of  heat  in  one  dimension  in  a  gas.  Here  the 
equation  is 

d2ip  2  dip  1  d2ip 
dx1  °  dt  +  c2  dt 2 

(see  Morse  and  Feshbach  [1953],  pp.  865-9).  In  this  equation,  ip  is  absolute  temperature,  a 
is  a  property  of  the  medium  and  c  is  the  velocity  of  sound  in  the  medium.  This  equation  may 
disturb  some  who  are  used  to  the  equation  for  transmission  of  heat  being  the  traditional 
Heat  Equation :  ipt  =  (1  /a2)ipxx.  The  usual  heat  equation  implies  that  heat  moves  infinitely 
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fast  (John  [1978],  ch.  7),  which  is  not  physical.  The  equation  given  here  has  the  property 
that  heat  can  move  no  faster  than  the  sound  speed  of  the  material  (as  is  easily  seen  by 
writing  the  equation  in  characteristic  form).  The  reason  the  heat  equation  is  such  a  good 
approximation  to  the  heat  transmission  problem  is  that  it  is  the  limit  as  c  -*  oo  of  the  correct 
equation.  That  is,  the  transmission  of  heat  is  so  slow  relative  to  c  that  the  finite  speed  of 
propagation  has  little  effect  on  the  solution. 

We  note  that  Morse  and  Feshbach  [1953]  provide  a  solution  to  this  equation  by  means 
of  Green’s  functions.  This  solution  is  not,  however,  very  enlightening  or  generally  useful. 

Neutron  transport  provides  another  source  for  these  equations.  The  actual  equations 
are  complicated  integral  equations;  see  Bell  and  Glasstone  [1970]  and  Richtmyer  and 
Morton  [1967],  chapter  9.  A  typical  approximation  to  these  equations  assumes  that  the 
neutrons  are  confined  to  a  small  discrete  set  of  energies.  In  this  case,  a  coupled  system  of 
equations  similiar  to  the  telegrapher’s  equation  arises.  The  lower  order  terms  in  this  case 
model  the  loss  of  neutrons  due  to  absorption. 


Historical  Background 

The  historical  background  for  this  thesis  consists  of  two  parts.  First  is  a  review 
of  the  work  on  lower  order  terms  in  constant  coefficient  hyperbolic  partial  differential 
equations  Second  is  an  outline  of  some  methods  for  studying  the  approximate  solution  to 
these  equations  analytically.  As  an  introduction,  we  discuss  the  two  main  effects  that  the 
presence  of  lower  order  terms  can  have  on  the  solution. 

The  first  effect  of  a  lower  order  term  is  the  obvious  one:  changing  the  equation.  In 
problems  from  physics,  the  lower  order  term  represents  a  form  of  dissipation.  For  example, 
consider  the  equation  ut  =  ux  4-  au.  It  is  easy  to  see  that  the  solution  to  this  is  u(x}t)  ~ 
u(x  +  t,  0)c“‘.  Thus,  if  a  <  0,  the  lower  order  term  damps  the  solution. 

Another  possible  effect  of  lower  order  terms  is  the  coupling  that  they  provide  among 
equations.  Consider  the  hyperbolic  system 


du  du 

-  =  4-+«„. 


(12) 


Here,  A  and  I )  are  constant  matrices,  and  u  is  the  vector  of  dependent  variables.  For  (1.2) 
to  be  hyberbolic,  the  eigenvalues  of  A  must  be  real.  This  is  the  general  constant  coefficient 
case,  since  by  the  assumption  of  hyperbolicity,  any  matrix  multiplying  the  <)u/dt  term  can 
be  inverted  and  moved  into  A  and  R.  Further,  we  assume  that  A  is  diagonalizable  with 
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real  eigenvalues  so  that  the  system  will  be  hyperbolic.  Now,  if  A  and  //  do  not  commute, 
so  that  they  are  not  simultaneously  diagonalizable,  this  will  cause  coupling,  (if  A  and  //  do 
commute,  the  coupling  is  apparent  but  not  real.  This  will  be  clear  in  the  following.)  To  see 
this,  make  a  change  of  variable  Tu  —  v,  where  TAT~l  =  I),  I)  a  diagonal  matrix.  Then 
we  are  left  with 


<)v 

dt 


D -f  TUT 

OX 


V. 


This  equation  is  as  simple  as  possible  while  still  exhibiting  coupling.  The  77/7'  1  term  is 
not  diagonal  if  A  and  //  do  not  commute  and  thus  couples  the  components  of  v. 

There  has  not  been  much  interest  in  the  effect  of  lower  order  terms  in  approximations 
to  constant  coefficient  hyperbolic  partial  differential  equations  until  recently.  A  major 
reason  for  this  is  that,  for  the  continuous  problem,  well-posedness  is  independent  of  the 
presence  of  lower  order  terms  (Thomee  [1969],  sec.  2).  This  is  also  true  for  difference 
approximations.  The  lower  order  terms  are  unimportant  in  the  question  of  stability  of  finite 
difference  schemes,  both  for  the  Cauchy  problem  (Thomee  [1969],  sec.  5)  and  for  the  initial 
boundary  value  problem  (Gustafsson,  Kreiss,  and  Sundstrom  [1972],  Thm.  4.3). 

There  has  been  some  work  on  equations  with  lower  order  terms.  Early  work  by 
Apelkrans  [1968]  on  scalar  equations  of  the  form 


du 

dt 


Oil 

=  pOM)^  +  <MM) 


u 


provides  bounds  on  the  size  of  errors  made  by  difference  approximations  at  a  step 
discontinuity.  These  bounds  are  sharp  for  the  general  problem,  but  may  be  somewhat 
pessimistic  for  specific  problems.  The  results  are  based  on  some  concepts  from  stability 
theory  for  difference  approximations  to  hyperbolic  problems  (see  Kreiss  and  Lundqvist 
[1968]),  and  are  difficult  to  extend  to  systems  of  coupled  equations.  Apelkrans’  results 
show  quite  generally  that  for  scalar  (or  uncoupled)  equations,  lower  order  terms  have  little 
effect  on  the  accuracy  of  difference  approximations.  Later  work  by  Brenner  and  Thomee 
[1971]  sharpened  these  results;  their  arguments  are  also  based  on  stability  theory  (Brenner 
and  Thomee  [1970]). 

Recent  work  by  Majda  and  Osher  [1977]  shows  that  lower  order  terms  can  have  an 
enormous  effect  on  the  error  in  difference  approximations  by  providing  the  kind  of  coupling 
discussed  above.  They  analyze  the  problem 
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dt\U2J  VO  -l/dzUa/ +  V-l  o)\uJ 

fl.  if  z  >  0 
v.\{x,0)  =  < 

lo.  if  X  <  0 


U2(z,  0)  =  0. 

Note  that  the  matrices  in  this  equation  do  not  commute.  They  show  that  if  U  is  a  dissipative 
finite  difference  approximation  to  (1 .3)  which  is  at  least  second  order  accurate  and  satisfies 
certain  technical  conditions  (which  are  satisfied  by  common  methods),  then  for  all  6  >  0 
there  exists  a  C6  >  0  such  that 

max  jC/  —  iz|  <  C(h 2 

x,teRf 

where  R&  =  {(z,f)  :  |z|  <  t  and  |z  ±  t\/t  >6}  and  0  <  6  <  t  <  7'0,  7’o  the  maximum 
time  for  the  problem.  This  result  is  the  best  possible  for  arbitrarily  accurate  dissipative 
schemes.  It  is  this  fact,  that  the  estimate  is  sharp,  which  is  the  interesting  feature.  This 
shows  that  over  a  large  region,  any  reasonable  dissipative  method  of  order  at  least  2  is 
only  second  order  accurate.  Note  that  the  solution  to  the  continuous  problem  is  smooth 
( C°° )  away  from  the  fronts  (at  x  —  —t  and  x  —  t)  since  the  initial  data  is  smooth.  The 
results  of  Majda  and  Osher  [1977]  do  not  hold  in  the  vicinity  of  the  step  discontinuity. 

We  now  discuss  some  of  the  methods  for  studying  the  behavior  of  difference  ap¬ 
proximations  to  hyperbolic  partial  differential  equations:  stability  theory.  Fourier  transforms 
of  the  difference  approximation,  model  equations,  and  an  ad  hoc  approach.  Methods  using 
stability  theory  are  very  powerful.  General  estimates  on  the  accuracy  of  the  solution  can  be 
found  for  a  wide  range  of  problems  with  this  method.  Examples  are  Lax  [1961]  for  Cauchy 
problems  and  Gustafsson  [1975]  for  initial-boundary  value  problems.  The  only  drawback  of 
this  method  is  that  it  requires  the  solution  to  be  sufficiently  smooth.  These  methods  do  not 
say  anything  about  the  behavior  of  the  difference  scheme  near  a  discontinuity.  The  results 
of  Apelkrans  [1968]  and  Brenner  and  Thomee  [1971]  mentioned  above  are  an  exception. 

Most  work  on  the  behavior  of  difference  approximations  to  hyperbolic  problems  by 
Fourier  transform  techniques  has  been  restricted  to  scalar  problems  (see  Ffedstrom  [1975] 
and  Apelkrans  [1968]).  This  is  in  contrast  to  stability  theory  for  difference  approximations, 
where  general  systems  of  linear  equations  can  be  treated  (see  Richtmyer  and  Morton 
[1967],  ch.  4  and  Coughran  [1980]).  The  advantage  of  the  Fourier  technique  is  that  it  can 
treat  nonsmooth  solutions  (such  as  step  discontinuities). 
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For  constant  coefficient  problems,  methods  using  Fourier  transforms  of  the  difference 
scheme  are  practical.  The  inverse  Fourier  transform  is  then  evaluated  by  the  method  of 
steepest  descent.  Hedstrom  and  Osterheld  [1980]  present  a  nice  example  of  this  procedure 
for  a  different  problem.  Pearson  [1969],  using  a  slightly  different  approach,  studies  the 
behavior  of  difference  schemes  for  ult  =  c2uxz. 

The  model  or  modified  equation  approach  is  another  way  to  study  the  behavior 
of  difference  approximations  to  partial  differential  equations.  In  this  method,  a  partial 
differential  equation  whose  solutions  model  those  of  a  difference  approximation  to  a  "'tnpler 
differential  equation  are  constructed.  The  model  equation  is  then  studied,  rather  than  toe 
original  difference  equation.  This  approach  is  discussed  by  Warming  and  Hyatt  [1974]  and 
in  a  review  by  Chin  and  Hedstrom  [unpublished].  As  we  will  use  this  approach  to  study  the 
behavior  of  general  difference  schemes  for  (1.1),  we  will  outline  it  below. 

The  form  of  the  model  equation  for  a  hyperbolic  problem  u<  +  l, a  —  0,  where  /  is  a 


differential  operator,  is 


On  _  .  p  dp+ 1  <t 

in  +Lu-Crh  0x p»T  +  ,»/i  ■’ 


where  p  is  odd  and  </  is  even.  Given  a  hyperbolic  differential  equation  ut  i-  Lu  —  I)  and  a 
difference  scheme,  expand  the  finite  differences  in  a  formal  Taylor  series  to  get 

du  °° 

+  Ij71  —  hrnl>m+  lu-  (f-5) 

m  — -  I 

Here,  /Jm  is  a  homogenous  differential  operator  in  <)x  and  dt  of  order  m.  Use  the  ansat/ 

=  V  cmhm  -  ■  — r  (1.6) 

Dt  ^  Oxm  +  ' 

m  —  t 

to  replace  all  derivatives  with  respect  to  t  in  (1.5)  with  x  derivatives,  and  then  match  up  the 
terms  in  (1.4)  and  (1.6),  dropping  higher  order  terms,  to  find  the  cm  in  (1.4).  This  gives  the 
model  equation. 

The  behavior  of  the  model  equation  can  be  studied  with  Fourier  transform  techniques. 
By  using  the  model  equation  approach,  we  can  avoid  some  of  the  algebra  required  when 
using  the  discrete  Fourier  transform.  More  important,  the  model  equation  approach  allows 
us  to  identify  the  key  features  of  the  approximation,  allowing  us  to  consider  general 
difference  approximations.  This  approach  has  been  taken  by  a  number  of  authors.  Using  the 
model  equation  formulation,  Hedstrom  [1975]  analyzed  the  behavior  of  step  discontinuites 


for  the  equation  ut  —  ux.  Work  hy  Ohm  [1975]  analyzes  the  detailed  behavior  near  the 
front  for  the  wave  equation.  Serdjukova  [1971]  considers  the  behavior  near  a  step  of  a 
general  scalar  difference  scheme  for  both  explicit  and  implicit  methods. 

Finally,  we  mention  an  approach  taken  by  Orszag  and  Jayne  [1974],  They  consider 
ut  +u,  =  0  and  look  for  solutions  with  continuous  derivatives  of  order  up  to  n  and 
discontinuous  derivative  of  order  n  +  I.  Through  a  clever  choice  of  initial  data  with  these 
properties,  the  difference  between  the  appoximation  and  the  true  solution  can  then  be 
estimated.  Chin  [1974]  shows  that  this  analysis  is  similiar  to  the  Fourier  analysis  above  and 
that  their  results  are  explained  by  model  equation  analysis. 


Brief  Introduction  to  Asymptotic  Estimation  of  Integrals 

We  present  here  a  brief  outline  of  the  method  of  steepest  descent,  on  which  most  of 
our  results  are  based.  More  complete  and  rigorous  descriptions  of  this  method  can  be 
found  in  Bleistein  and  Handelsman  [1975],  Erdelyi  [1956],  and  Olver  [1974],  trdelyi  [1956] 
in  particular  is  a  good  first  introduction  to  this  material.  We  will  first  introduce  the  method 
of  stationary  phase  because  it  may  be  more  familiar. 

We  will  need  to  find  approximations  to  integrals  of  the  form 

•  f  oo 


/ 


!/U! 


dx 


(1.7) 


when  t  is  large.  If  <j>[x )  =  i>J>(x),  and  ij>(x)  is  real  then  we  can  use  the  method  of  stationary 
phase.  The  reason  for  this  name  will  become  clear  in  what  follows. 

We  can  think  of  i!j(z)  as  a  frequency  of  oscillation  of  the  integrand,  if  we  assume  that 
<j(x)  is  slowly  varying.  Now,  look  for  the  extrema  of  y(x)  (these  are  points  of  stationary 
phase).  Assume  there  is  only  one  extremum,  x0.  Expand  both  rj[x)  and  >l'(x)  around  x0: 


g{x)  —g[x0)  +  (x  -  z0)^ (x„)  +  •  •  • 
rj>(x)  =iP(x „)  4-  (x  -  X(,)2^’(x0)  + 


(1.8) 


Using  these  expansions,  we  can  rewrite  (1.7)  as 

i*  +  oo 


Kz,,),^o)  J  oxp((f(x  -  x„)2V’"(x,)))</x  +  •  •  j. 

The  integral  in  the  formula  above  can  be  evaluated  exactly. 

This  works  because,  away  from  xo,  the  integrand  oscillates  rapidly  (since  t  is  large) 
and  hence  nearly  cancels  out.  Only  near  .r()  does  the  integrand  oscillate  slowly  in  x.  We 
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can  get  more  terms  in  the  expansion  simply  by  keeping  more  terms  in  the  expansions  (1.8) 
of  g  and  il>.  There  are,  however,  a  number  of  restrictions  to  remember.  First,  it  is  very 
important  that  g  vary  slowly  on  the  scale  of  exp (iip(x)t).  For  example,  if  g  has  a  pole  near 
an  extremum  of  ip(x),  the  expansion  (1.8)  will  be  inaccurate.  Also,  in  the  case  where  ip(x) 
has  several  extrema,  they  must  be  well  separated.  Otherwise,  there  may  be  an  interval, 
rather  than  a  point,  where  the  integrand  oscillates  slowly. 

What  if  4>(x)  or  xo  is  not  real?  Then  we  must  go  to  the  method  of  steepest  descent. 
We  extend  <A(x)  (we  will  not  use  i/'(x)  again)  to  be  a  complex  valued  function  <p(z)  over  the 
complex  plane.  This  is  usually  a  trivial  step  in  practice.  Since  <t>(z)  is  now  complex  valued, 
it  is  not  possible  to  say  that  oscillates  at  all  in  general.  Further,  if  4>(z)  is  analytic,  it 
has  no  minima  or  maxima  (by  the  maximum  modulus  principle). 

We  consider  the  following.  If  we  could  find  a  path  of  integration  I'  such  that 
|e*|/max,er|e*|  >  c  only  along  a  few  short  intervals,  then,  just  as  before,  we  could 
expand  g  and  <t>  about  points  on  those  intervals.  Let  us  assume  that  there  is  only  one  such 
point,  and  call  it  z0.  We  can  then  use  the  Cauchy  integral  theorem  to  deform  the  original 
path  of  integration  to  the  new  path  l\  Then  as  t  increases,  the  contribution  to  the  integral 
from  points  away  from  z0  will  decay  exponentially. 

A  good  way  to  look  at  this  is  in  terms  of  a  topographical  map.  Let  x  +  iy  =  z  be  the 
map  coordinates,  and  let  9?(<£(z))  =  real  part  of  <j>{z)  be  the  height  of  the  “ground.”  This 
map  will  normally  have  lots  of  valleys,  hills,  and,  most  importantly,  saddle  points  (maybe 
only  one).  A  simple  saddle  point  is  like  a  mountain  pass  —  at  the  saddle  point,  you  can 
either  go  downhill  (along  the  pass)  or  uphill  (at  right  angles  to  the  pass).  In  our  application, 
the  path  typically  starts  in  a  valley  at  x  =  —  oo,  proceeds  over  one  or  more  mountain 
ranges,  and  ends  in  a  valley  at  x  =  oo.  When  going  over  the  mountain  range,  the  best 
path  to  take  is  one  going  through  the  saddle  points.  The  path  then  stays  low  except  for 
a  few  peaks  where  it  rises  quickly  to  go  though  a  pass  (saddle  point)  and  then  falls  again 
quickly.  This  is  what  we  want,  since  the  exponential  will  have  small  magnitude  everywhere 
except  in  a  small  neighborhood  of  each  pass  (saddle  point).  This  explains  why  the  method 
is  called  the  method  of  steepest  descent,  as  the  path  descends  as  quickly  as  possible  from 
each  pass. 

Given  this  geometric  viewpoint,  it  is  easy  to  find  conditions  for  the  passes  and  for  the 
direction  of  the  path.  These  conditions  are  derived  in  Bleistein  and  Handelsman  [1975], 
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chapter  7.  Basically,  an  n-th  order  saddle  point  zq  is  a  point  where  the  first  n  derivatives  of 
4>(z)  are  zero  and  the  n  +  1st  derivative  is  non-zero.  This  condition  is  the  same  as  for  the 
method  of  stationary  phase  discussed  above.  Then  it  can  be  shown  (see,  e.g.,  Bleistein  and 
Handlesman  [1975],  Thm  7.1)  that  the  direction  of  the  path  of  steepest  descent  from  an  n-th 
order  saddle  point  at  the  saddle  point  is  given  by  the  following  rule:  If  dn+l<j>{zo)/dzn+l  — 
ae'a,  then  the  directions  of  steepest  descent  are  given  by  0  =  (— a  +  (2 p  +  l)7r )/n,  for 
p  =  0,l,...,n-  1. 

The  above  sketch  of  the  method  of  steepest  descent  ignores  a  number  of  possible 
difficulties.  Most  important  for  our  purposes  is  what  happens  if  g  is  not  slowly  varying 
or  if  several  saddle  points  come  together.  The  proper  way  to  handle  this  is  discussed  in 
section  2;  basically  it  amounts  to  using  better  expansions  for  g  and  <j>.  The  difficulty  with 
these  more  accurate  expansions  is  that  the  integrals  in  the  expansion  may  no  longer  be 
recognizable  as  known  special  functions.  Instead,  they  may  define  new  special  functions 
whose  behavior  must  be  investigated.  We  will  encounter  a  new  class  of  special  functions 
in  sections  3  and  4. 
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2.  The  Exact  Solution 


In  this  section  we  find  the  exact  solution  to  the  model  problem 


ut  =  - v , 


vt  =  —  Ux  —  V 


u(x,  0)  =  v(x,  0)  = 


1, 

I 

2’ 

0, 


if  x  <  0 
if  x  =  0 
if  x  >  0. 


(2.1) 


We  are  only  interested  in  u  near  the  front  (x  =  t).  Our  solution  will  be  valid  only  in  this 
region. 

Our  purposes  for  finding  the  exact  solution  to  the  model  problem  are  manifold.  First, 
we  will  need  the  exact  solution  to  compare  with  the  solutions  of  approximations  to  the 
problem.  Second,  in  finding  an  asymptotic  expression  for  the  exact  solution  we  will 
demonstrate  many  of  the  techniques  that  we  will  use  in  later  sections.  Most  important,  we 
will  use  the  analysis  in  this  section  to  motivate  definitions  and  analysis  in  later  sections. 

We  point  out  that  for  x  >  t,  u(x,  t)  =  v(x,  t)  =  0  from  the  theory  of  characteristics 
(John  [1978]).  Unless  otherwise  stated,  all  limits  x/t  — ►  1  are  from  the  left  in  the  x-t  plane. 

The  analysis  will  proceed  as  follows.  First,  we  Fourier  transform  (2.1)  in  space  to  get 
a  system  of  two  ordinary  differential  equations.  As  these  are  linear  constant  coefficient 
equations,  it  is  easy  to  solve  for  the  inverse  transform  of  the  solution.  The  inverse  transform 
does  not  represent  any  known  special  function,  so  we  will  use  techniques  from  saddle  point 
analysis  to  represent  the  inverse  transform  near  the  front  x  =  t  as  an  asymptotic  series  in 
terms  of  modified  Bessel  functions. 


Integral  Form  of  the  Solution 
The  Fourier  transform  of  (2.1)  is 

ut  =  —  t£v 
vt  =  -i£u  —  v 

where  £  is  the  dual  variable.  In  matrix  form,  this  is 


dt 


(22) 
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The  eigenvalues  X  of  the  matrix  are  given  by  \(\  +  1)  +  £a  =  0,  or 


X±U)  =  -i±iv/4^TT. 


The  branch  cut  for  the  square  root  is  the  interval  in  the  (  plane.  This  is  the 

standard  branch.  The  unnormalized  eigenvectors  of  the  matrix  are 


Therefore,  the  solution  to  (2.2)  is 

where  a±(£)  are  determined  by  the  Fourier  transform  of  the  initial  conditions. 

Now  that  we  have  the  solution  to  the  transform  of  (2.1),  we  can  compute  u  and  v  by 
inverting  the  transform.  The  solution  for  u  is  thus 


^ t] = h  /:;  mw + x+)cx+t + a-(f ){i + x-)ex"‘]e*x<  <2-3> 


Here  we  choose  the  path  of  integration  to  pass  over  the  branch  cut  between  -1/2  and  1/2 
(an  arbitrary  choice). 

To  determine  ot  we  apply  the  initial  conditions  to  the  solution  of  the  differential 
equation.  At  t  =  0,  we  have 

o+tf)(l  +  X+)  +  +  X_)  =  f{() 

-m+(0e-*a_(fl*  ==/(*) 

where  u(£,  0)  =  v(£,  0)  =  /(£).  These  imply 


“-(f,=dfe(,-l(i+x+)) 


a+(C)=*Ai)-a_(0- 

For  the  initial  data  given  by  (2.1),  we  have  /(()  =  ft  +  f  6(0- 


(2.4) 


Solution  near  the  Front 

We  have  now  completely  specified  the  solution  to  (2.1).  We  need  only  evaluate  the 
integral  (2.3)  in  order  to  determine  the  solution.  However,  there  is  no  Known  special 


function  that  this  integral  represents.  To  uncover  the  solution,  we  will  turn  to  asymptotic 
evaluation  of  the  integral  by  the  method  of  steepest  descent  outlined  in  the  introduction. 
We  will  outline  the  steps  again,  as  we  will  see  that  there  are  new  complications:  converging 
saddle  points  and  a  pole  in  g.  We  write  (2.3)  as  the  sum  of  two  integrals  of  the  form 

I  gU)e+{()tdt. 

Here,  g  does  not  contain  an  exponential  factor  in  t  or  x. 

The  analysis  of  (2.3)  takes  place  in  five  steps.  First,  the  saddle  points  of  <t>  are  found. 
Next,  the  path  of  integration  is  deformed  to  the  steepest  descent  path  through  some  of  the 
saddle  points.  We  will  call  these  the  important  saddles.  This  path  will  not  necessarily  pass 
through  all  of  the  saddle  points.  Care  must  be  exercised  at  this  point  if  the  original  path  of 
integration  passes  through  any  singularities  of  the  integrand,  if  the  original  path  does  pass 
through  some  simple  poles  and  the  new  path  does  not,  half  the  residues  at  those  points 
will  have  to  be  added  to  the  integral.  Third,  a  change  of  variable  is  made  to  map  <j>  onto 
a  simpler  function  0  that  has  the  same  number  of  important  saddle  points  as  <j>  has.  Next, 
g  is  expanded  in  a  power  series  about  the  saddle  points  of  0  and  any  critical  points  of  g 
(in  our  case,  we  expand  only  about  the  critical  points  of  g;  this  will  be  discussed  later). 
Finally,  the  resulting  integral  is  evaluated.  In  this  section,  these  integrals  can  be  evaluated 
in  terms  of  known  special  functions.  If  this  is  not  the  case,  as  in  the  next  section,  then  the 
integrals  represent  a  new  special  function  which  must  be  examined.  The  result  we  now 
prove  is: 

Theorem  1:  For  x  <  t  and  x/t  «  1,  the  solution  for  u  in  (2.1)  has  the  asymptotic 
representation 

“(*»0  =  £  =  \e~ 1,2  £  ([~)  M^\/l-w2).  (2.5) 

i/=0  i/=0  '  ' 

Proof: 

We  first  find  the  saddle  points.  Let  w  =  x/t,  and  define  0±(£)  =  X±(£)  +  tw£.  Then 
(2.3)  can  be  written  as 

J g+et+'di  +  J g.et-'dt. 

Each  of  these  is  in  the  form  discussed  above.  The  condition  for  the  saddle  points  of  <f>±  is 

%(0  =  ±  — ^ =  0.  (2.6) 
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The  saddle  points  are 


e*  =  ±- 


w 


„ _  (2.7) 

2\ZwJ  -  1 

It  is  easy  to  show  that  </2^±/d^a  0  for  u  ^  1,  so  these  saddles  are  simple  for  w  ^  I. 

At  (J  =  1 ,  the  saddle  points  are  both  at  infinity,  and  are  of  infinite  order. 

Now,  we  must  be  careful  here  as  the  formula  (2.7)  is  the  result  of  solving  a  quadratic 
formed  by  squaring  (2.6).  This  means  that  we  must  check  to  see  which  choice  of  sign,  if 
any,  is  a  solution  to  the  equation  (2.6).  This  is  equivalent  to  determining  on  which  side  of 
the  branch  cut  of  the  square  root  the  solutions  (2.7)  lie.  Note  that  there  are  four  conditions 
to  check  —  in  both  <f>~  and  4>+. 

By  substituting  (2.7)  into  (2.6),  we  can  show  that  for  w  >  0,  <t>+  has  no  saddle  points 
and  4>-  has  two  saddle  points.  For  w  <  0  the  situation  is  reversed. 

It  is  thus  natural  to  break  the  integral  (2.3)  into  two  parts 

U+  =  ‘Ll-00  "  +  +  X+)e*+‘d£ 


-oo 
r  +  oo 


i  f+°° 


(2.8) 


dt. 


We  show  in  Appendix  A  that  u+  =  0  for  w  >  0.  We  therefore  concentrate  on  the  behavior 
of  u_,  in  particular  for  w  «s  1. 

From  (2.7),  we  see  that  as  w  -*  1,  the  saddle  points  go  to  ±ioo.  The  major  contribution 
to  the  integral  will  come  from  two  places.  One  is  the  region  of  the  saddle  points;  the  other 
is  the  pole  in  the  integrand  at  £  =  0.  The  path  may  be  deformed  away  from  the  (£  ±i)-,/2 
singularities  in  a_(l  +  X_)  without  changing  the  value  of  the  integral,  as  these  singularities 
are  weaker  than  poles. 

We  first  deform  the  path  away  from  the  pole  at  the  origin.  To  correct  for  this  (since  the 
integral  is  the  Cauchy  Principle  Value)  we  add  7tj  times  the  residue  at  £  =  0.  From  (2.4) 
and  (2.8)  it  is  easy  to  see  that  the  residue  is  i/4n  (for  the  path  of  integration  passing  over, 
rather  than  under,  the  branch  cut),  yielding  a  term  of  -1/4  to  be  added  to  the  integral  over 
the  new  path. 

We  can  now  consider  the  effects  of  the  saddle  points.  Since  there  are  two  saddle 
points,  we  will  need  to  consider  a  uniform  approximation  to  the  integral  which  models 
the  confluence  of  two  saddle  points  to  a  saddle  point  of  infinite  order.  Such  a  case  is 
considered  by  Bleistein  and  Handlesman  [1975]  in  section  9.5.  We  will  use  essentially  the 
same  approach  as  they  use  for  this  problem. 
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To  make  the  analysis  easier,  we  first  map  (  onto  1/z,  and  work  in  the  z-plane.  In  the 
z-plane,  the  saddle  points  are  converging  to  the  origin  as  w  -*  1. 

Expanding  <p_(l/z)  about  z  =  0,  we  see  that 

,  It  f.  ?  iu 

=  t  +  t 

Using  this  as  motivation,  we  define  a  change  of  variables  p  =  p(z)  by 

V»(p)  —iP  +  b+~  =  *-(7)-  (2-9> 

o  p  z 

Here,  a  and  b  are  numbers  independent  of  p  but  perhaps  depending  on  w.  It  is  important  to 
note  that  this  form  is  the  simplest  form  which  has  the  requisite  properties.  The  requirements 
that  ip  have  exactly  two  saddle  points  which  come  together  to  a  single  saddle  point  of  infinite 
order  completely  determines  the  simplest  form  of  ip. 

To  find  a  and  b,  let  z±  =  l/(±,  the  saddle  points  of  <p(£),  and  let  p±  =  p(z±)  be 
the  saddle  points  in  the  p-plane.  Since  p(z)  is  locally  1-1  (see  Guillemin  and  Sternberg 
(1977],  page  441-2),  we  have 

d<P_ 

dz  dz  /  dp  t  o' 

8  ~  P* 

The  fact  that  p(z)  is  locally  one-to-one  everywhere  means  that,  in  particular,  it  is  one-to-one 
at  the  saddle  points.  Thus  we  must  have  dz/dp  ^  0  at  the  saddle  points.  Since  by 
definition  dp_/dz  =  0  at  the  saddle  points,  we  must  have 

*  oa 

8  p(*J8 

or  p2(z±)  =  — 8ia2.  Applying  these  two  equations  to  (2.9)  gives  ±2\/ij&a  +  b  =  <P-(£±)  or 

6  _*_(*+)+< My 

v»  ' 

Being  careful  with  the  branch  of  the  square  root,  it  is  easy  to  show  that  <P-((±)  =  -  J  ± 
i-v/w1  —  I,  and  hence  that 


Thus  V>(p).  the  canonical  form,  is 


V'(p)  = 


(2.10) 


Note  that  for  u  s=»  1,  this  is  very  close  to  the  z  -»  0  limit  of  <f>_. 

The  unexponentiated  part  of  the  integral  has  three  components.  First  is  the  /3(1  +  X_) 
from  (2.3).  Second  is  a  -1  /z2  term  coming  from  the  change  of  variables  £  =  1/z.  Finally, 
there  is  the  dz/dp  term  coming  from  the  change  of  variables  p  =  p(z).  To  determine  these 
in  terms  of  p,  we  need  to  find  a  representation  for  p  =  p(z).  We  also  need  to  find  dz/dp. 
In  a  general  case,  we  can  find  this  information  from  (2.9);  however,  in  this  case,  we  can 
find  the  map  explicitly. 

Now,  (2.9)  may  be  thought  of  as  a  quadratic  equation  in  z,  whose  coefficients  depend 
on  p: 

ip  «(^2  -  1)  _  *w  _  _»  A  _  f2 

8  2p  z  zV  4 


or 


'P  ,  («a-0  , 

.8  +  2  p 


)2 


Z* 

T 


With  some  algebraic  manipulation,  we  may  rewrite  this  as 


_  8p(u;  T  1) 

p2  +  4(w  =F  l)2 

To  pick  between  the  two  solutions  for  z,  we  insist  that  for  w  near  1  ,z«p,  with  both  p  and 
z  near  0.  This  requirement  just  says  that  p  =  p(z)  is  one-to-one  in  the  neighborhood  of 
the  saddle  points  near  z  =  0.  This  is  satisfied  by  the  +  choice  in  T,  and  we  have  finally 


z  = 


2  P 


uj  +  1 


1 


1  + 


( 


— ? i .  y 

2(w  +  i)y 


It  will  be  useful  to  expand  this  about  p  =  0;  the  result  is: 


2p 

u  +  1 


P 

2(w  +  1) 


(2.11) 


Differentiating  (2.11),  we  find  that 


dz 

dp 


2 

w  +  1 


OO 


Ei-'ni 


+ 


P 

2(w  +  1) 


2  3p2 


>ir 


I 


This  is  the  appropriate  form  of  dz/dp,  since  the  saddle  points  are  converging  to  the  origin, 
and  the  unexponentiated  term  has  a  pole  at  p  =  0.  (We  have  used  the  form  recommended 
by  van  der  Waerden  [1951]  for  saddle  points  near  a  pole.  Other  expansions  are  possible; 
see  Bleistein  [1966]  and  Olver  [1974].) 

Next,  we  find  the  unexponentiated  part  of  the  integrand.  It  is  convenient  to  consider 
<*_(1  +  X_)(— 1/z2)  =  h(p)  together.  Neglecting  the  6-function  term  in  a_(l  +  X+),  this  is 

-  v=('  - 

Now,  we  need  to  convert  this  to  the  p  variable.  We  handle  the  square  root  first. 


16p2(w  +  l)2 
(p2  +  4(w  +  l)2)2 


4((V  -)-  l)2  —  p2 
4  (ui  +  1)2  +  p2 


where  the  bt  «.t>ch  ol  'she  square  root  was  chosen  for  p  small.  Using  this  formula,  we  find 
after  some  manipulation  that 


,  _  p2  +  4(w  +  l)2  (  2 i(w  +  1)\ 

‘w -  ?){'  - ~ir-} 

Multiplying  this  by  the  last  part  of  the  unexponentiated  term,  dz/dp,  we  get 

i  (i  _  2t-(«+in  i 
4(w  +  l)V  P  /  p 2 


4(w+  1)2 


+  1 


Finally,  as  a  power  series,  this  is 


4( 


_1 _ y  (  -V 


(2.12) 


We  are  now  ready  to  write  u  as  an  asymptotic  series.  Combining  (2.10)  and  (2.12),  we 

have 


u(x,  t)  —  - 


87r(w  +  1) 


Again,  we  have  deformed  the  path  of  integration.  The  6(f)  term  in  /  which  we  have  been 
neglecting  cancels  the  —  |  term  from  the  pole  at  £  =  0  (which  is  z  =  oo).  This  pole  is  not 
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present  in  our  expansion  of  the  unexponentiated  term  because  we  broke  the  integral  u_  in 
(2.8)  into  two  parts:  one  part  about  the  saddle  points  near  z  =  0.  and  the  other  part  near 
the  pole  at  a  =  oo.  We  can  deform  the  path  away  from  the  pole  at  z  =  oo,  including  the 
contribution  from  the  pole.  This  leaves  us  with  just  the  integral  above.  (Note  that  there  is 
also  a  pole  at  p  =  0.  This  pole  is  handled  as  part  of  the  integral  —  see  the  second  integral 
identity  in  Appendix  A.) 

The  first  term  in  this  expansion  is 


*0 


dp. 


The  second  integral  identity  in  Appendix  A  tells  us  that  this  is 


“o  =  ^r‘/2|o(^Vl  -w2), 


where  l„  is  the  regular  special  Bessel  function  of  order  u.  This  is  a  good  approximation 
for  w  less  than  and  near  1. 

The  next  term  is  just 


From  the  formula  in  Appendix  A,  this  is 


This  process  may  be  carried  out  repeatedly;  this  completes  the  proof.  | 

It  is  interesting  to  consider  the  behavior  of  the  first  few  terms  in  the  approximation 
(2.5)  near  the  front  (w  s=»  1).  First  we  note  that  for  w  =  1,  the  v  >  0  terms  are  all  zero 
since  l»,(0)  =  0  for  u  >  1  (see  Abramowitz  and  Stegun  [19721).  ,n  general,  at  w  =  1,  the 
u  >  k  terms  are  0  for  u{z,  <)  and  its  first  k  derivatives  in  z.  Thus,  the  first  k  terms  give  u 
and  its  first  k  -  1  derivatives  in  x  exactly  at  t  =  x.  We  expect,  therefore,  that  the  first  few 
terms  would  give  us  a  good  approximation  to  the  solution  for  x  *=w  t. 


Asymptotic  Behavior  near  the  Front 

We  can  get  a  better  idea  of  behavior  of  the  solution  near  the  front  x  =  t  by  looking 
at  the  small  argument  behavior  of  \u{z).  Using  the  small  argument  form,  (2.5)  gives 
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«(*,<)  =  ^-',,(i  +  lt£  +  ‘ 


2_22 


(t  ~  x)2 
+  —  + 


16  32 

£>,-*/»(,+ I +  ...). 

tfz  8  v  2  / 

It  is  clear  from  the  power  series  for  1^  that  uu  =  6*,(f  -  x)v  cxp(-</2)  +  0((f  - 
x)l/+2)  exp(— </2)  for  u  w  1.  At  t  —  1  and  x  =  1  these  give  u  =  £e_1/2  and  w*  = 
exactly.  We  can  compare  these  with  the  computed  solution  given  in  figure  1. 
The  method  used  for  computing  the  exact  solution  is  given  in  Appendix  B. 


Figure  1 :  Exact  solution  for  u  for  t  =  1. 

As  we  can  see  from  figure  1,  these  match  our  results.  In  particular,  note  that  the  height 
of  the  front  is  roughly  0.3,  which  agrees  with  the  computed  value  of  «=»  0.303. 

Also,  the  slope  near  the  front  can  be  estimated  as  0.3  —  0.41  =  —0.11,  which  matches  the 
calculated  value  of  — 3e—  1 /2/16  «  —01 14. 

Finally,  we  discuss  the  range  of  validity  of  the  expansion  (2.5).  The  limitation  comes 
from  the  fact  that  the  expansion  (2.12)  is  valid  only  for  \z\  <  2.  Thus,  as  w  decreases  from 
1,  we  can  expect  the  ( z  ±  2)~ 1|/2  terms  to  have  an  increasing  effect.  As  a  rule  of  thumb, 
since  |*J  =  2  at  w  =  (l/2),/l2  «  0.707,  the  approximation  will  be  valid  for  a  range  of  u 
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near  1,  such  as  1/2  +  l/2(l/2)1^  <  <*/  <  l.  If  we  wanted  an  expansion  that  was  valid 
over  a  greater  range  of  w,  we  would  have  to  either  use  an  expansion  of  «_(1  +  X_)  with 
a  greater  radius  of  convergence  or  use  the  method  of  matched  asymptotics.  In  this  case, 
the  method  of  matched  asymptotics  would  be  the  more  fruitful  approach.  In  this  approach, 
as  the  saddle  points  approach  the  branch  cut,  the  path  of  steepest  descent  begins  to  be 
affected  by  the  branch  cut.  At  the  same  time,  the  saddle  points  are  moving  apart  and 
hence  interacting  less  and  less  strongly.  Thus  we  could  use  an  expansion  valid  for  two 
saddle  points  approaching  a  branch  cut.  In  the  transition  region,  we  would  match  these 
two  asymptotic  expansions.  Since  we  are  interested  in  the  solution  only  near  the  front 
w  =  1 ,  we  will  not  consider  either  of  these  further. 
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3.  Behavior  of  Difference  Approximations  near  the  Discontinuity: 
Second  Order  Method  of  Lines 


in  this  section  we  analyze  a  particular  difference  approximation  to  the  model  problem. 
For  simplicity  in  the  discussion,  we  consider  the  method  of  lines  (or  semi  discrete)  second 
order  centered  difference  approximation 


dU, 

dt 


dt 


Vj+i-Vj-x 

2h 

2  h 


1 

2  * 


<Mo)  =v,-(o)  =  { 


if  j  <  o 
if  j  =  0 
if  j  >  0. 


(3.1) 


This  is  a  second  order  non  dissipative  difference  approximation  to  our  model  problem. 
The  space  step  size  is  h.  The  grid  is  not  staggered.  For  simplicity  of  analysis,  the 
undifferentiated  term  is  not  averaged.  Also  for  simplicity  of  analysis,  the  initial  conditions 
are  slightly  different  from  those  in  the  previous  section.  However,  this  will  not  affect  the 
behavior  of  the  difference  scheme,  since  both  the  equation  and  the  difference  scheme  are 
linear.  Of  course,  this  change  in  initial  data  changes  the  solution  everywhere  by  j  for  U 
and  by  ^ e~*  for  V  for  both  the  continuous  and  the  semi-discrete  problem,  but  we  can 
ignore  this  since  the  change  is  the  same  for  all  z. 

Our  analysis  will  illustrate  the  behavior  of  the  approximate  solution  to  U}(t)  near  the 
front  x  —  t.  We  are  particularly  interested  in  the  asymptotic  behavior  of  the  width  of  the 
front  as  h  -*  0.  In  addition,  we  will  show  that  the  behavior  of  the  difference  scheme  near 
the  front  is  very  similiar  to  that  of  the  wave  equation. 

Since  this  section  is  the  longest,  we  will  provide  a  more  detailed  summary  of  it.  We 
first  find  the  inverse  Fourier  transform  which  represents  the  solution  to  (3.1).  We  will  see 
that  the  inverse  transform  is  more  complicated  than  that  in  section  2;  in  fact,  its  canonical 
form  represents  a  new  special  function.  At  this  point  the  analysis  diverges  from  that  in 
the  previous  section.  We  take  up  the  question  of  the  width  of  the  front  in  (3.1)  (since  the 
discontinuity  in  the  solution  to  (2.1)  is  smeared  out  by  the  approximation  (3.1)).  First  we 
discuss  the  behavior  of  the  exponential  in  the  inverse  transform  integral.  This  leads  us  to 


a  definition  for  the  width  of  the  front.  We  then  go  on  to  state  and  prove  a  theorem  that  the 
width  of  the  front  is  0(h2/3).  After  the  proof  is  complete,  we  turn  to  the  canonical  form  of 
the  inverse  tranform.  From  the  canonical  form,  we  can  show  the  behavior  of  the  solution 
to  (3.1)  near  the  front  (mostly  through  graphs  displayed  in  section  5).  We  close  with  a  few 
secondary  issues  which  may  be  skipped  on  a  first  reading.  The  first  is  the  behavior  of  the 
solution  away  from  the  front.  The  other  discusses  the  size  of  a  term  that  was  neglected  in 
the  analysis. 


Integral  Form  of  Solution 

We  start  with  the  discrete  Fourier  transform  of  (3.1): 

i  — 

where  £  is  the  dual  variable.  Let  0  =  h£ .  In  matrix  form,  this  is 

dt\ V)  V— i"0  -*  A  VJ 

The  eigenvalues  X  of  the  matrix  are  given  by  X(X  +  I)  +  sin2  0/h2  =  0,  or 


(3.2) 


MO 


1  i 

2  *  2  ^ 


1  sin2  0 


-  1. 


The  branch  cut  for  the  square  root  is  the  interval  [-Arcsin(/i/2)//i,  Arcsin(/t/2)//t]  in  the  £ 
plane.  This  is  the  standard  branch.  The  unnormalized  eigenvectors  of  the  matrix  are 


Therefore,  the  solution  to  (3.2)  is 


(^)  -  Sin+0)<M  +  Mo(_<  8in'tf)eX-‘ 


where  MO  are  determined  by  the  discrete  Fourier  transform  of  the  initial  data.  The 
solution  for  U  is  then  just  the  inverse  transform 

*  +  ir/h 
-w  fh 


/•  +  »/* 

£/(M)  =  5-  /  MOO  +  x+)«x+*  +  MOO  +  x_) 

2ir  J-„/h 


(3.3) 


As  in  section  2,  the  path  of  integration  passes  above  the  branch  cut.  Let  4>±  =  X±  +  i(x/t). 
We  will  write  the  two  terms  in  (3.3)  corresponding  to  <f>+  and  <t> .  as  U+  and  V-  respectively, 
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as  in  (2.8).  We  note  that  this  formula  defines  U  for  all  real  x,  even  though  U  was  originally 
defined  only  for  x  =  jh  for  integer  j.  The  inverse  transform  (3.3)  provides,  however,  a 
natural  interpolation  to  real  x,  and  we  will  use  this  in  the  rest  of  this  section. 

To  determine  a±  we  apply  the  initial  conditions  to  the  solution  of  the  differential 
equation  (3.2).  At  t  =  0,  we  have 

«+(0(i  +  V)  +  «-U)(i  +  M=/(0 

,„.sin0  ,..sin0  iit\ 

a+(0-r-  +  Q-(£)~r-  =1^ 


where  t>(C,0)  =  V'(^O)  =  /(£).  These  imply 
<*-(£)  = 


-JL — 4*  (,  +  x,)") 

A  Sln0  J 

Jl*  1 


V  * 


or 


o.(0(*  +  M 


(3.4) 


<*♦100 +  M  =  -  , — — — 

2/^-1 

For  the  initial  data  given  by  (3.1),  we  have  /(£)  =  i7icot(0/2)/2. 

We  need  now  to  investigate  the  behavior  of  (3.3)  with  (3.4)  in  the  vicinity  of  x  =  t.  We 
start  by  defining  u>  =  x/t  and  0±(C)  =  X±(C)  +  iwf .  (Where  it  is  more  convienent,  we  will 
write  0(0 )  for  0(C).)  The  condition  for  the  saddle  points  is  then 


d<t>±.,.  2t  sin  0cos0  .  _ 

~fU)  =  ± - + =  0. 


(3.5) 


.  / i  sin2  6 


1 


The  solution  for  the  saddle  points  0±  is 


sin2  0±  — 


1  —  w 


2 


2 


(3.6) 


There  are  clearly  eight  solutions  C  to  this  equation  in  the  rectangle  (~n/h,  itjh )  X  (-00, 00). 
Let  0±  =  sin-1  »±,  with  0  <  0+  <  tt/2  and  3?0_  =  0,  0_/t  >  0.  (Slz  is  the  real  part  of 
z.)  Then  the  other  six  roots  are  —  0+,  it  —  0+,  -it  -F  0+,  — 0_,  it  +  0_,  and  7r  —  0_.  To  see 
which  of  these  solutions  are  zeros  of  0'_(C)  and  which  are  zeros  of  <f>'+((),  we  plug  (3.6) 
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into  (3.5).  Again,  we  must  be  careful  with  the  branch  of  the  square  root.  We  first  note  that 
sign(9J sin  0)  ==  sign(!R('l sin2 6/h 2  —  l)1/2).  Thus,  by  considering  the  sign  of  5?cos0,  we  see 
that  ^'_(0  has  roots  h(  =  ± 0 +  and  ±0_,  and  has  roots  ±7r  0 +  and  7r  ±  0_. 

Note  that  the  location  of  the  saddle  points  is  substantially  more  complicated  than  for 
the  integral  in  section  2.  This  complication  makes  it  impossible  to  obtain  an  asymptotic 
expansion  for  (3.3)  in  terms  of  known  special  functions.  We  will  have  to  be  content  with 
some  estimates  of  its  behavior.  Note  also  that  the  integral  of  the  <f>+  term  is  not  zero,  in 
contrast  to  the  continuous  case.  However,  we  will  show  later  that  it  is  of  the  same  size 
as  the  first  neglected  term  in  the  asymptotic  expansion.  For  the  moment,  though,  we  will 
ignore  the  </>+  term. 

Width  of  the  Front 

Perhaps  the  most  interesting  information  that  we  can  look  for  is  the  width  of  the  front. 
Serdjukova  [1971]  and  Hedstrom  [1975]  have  shown  that  for  the  scalar  wave  equation,  the 
width  of  the  approximation  to  a  step  discontinuity  by  a  difference  method  of  order  p  is 
0(/tpAp+,)).  yhis  smearing  out  of  the  front  can  be  thought  of  as  the  result  of  numerical 
dispersion  or  dissipation  or  both,  introduced  by  the  approximation.  In  this  section  we  will 
show  that  the  same  relation  holds  for  the  second  order  method  (3.1)  under  consideration. 
We  show  this  as  follows.  First,  we  study  the  path  of  steepest  descent.  This  will  show  us  that 
the  saddle  point  0_  controls  the  width  of  the  front.  Next,  we  use  perturbation  techniques 
to  find  0 _  and  <j>-(0-)  ’•  r  x/t  =  u  —  1  •+-  ch 7.  By  studying  the  behavior  of  <£_(0_)  as  a 
function  of  7,  we  will  see  that  the  width  of  the  front  is  indeed  0(h2/3). 

In  passing,  we  note  that  the  width  of  the  front  is  not  a  well  defined  concept.  This  is 
because  there  are  a  number  of  other  phenonema  which  interfere  with  any  measurments 
of  a  width  of  the  front.  For  example,  the  oscillations  which  are  present  in  almost  any 
difference  approximation  can  mask  the  width  of  the  front.  Despite  this,  it  is  sometimes 
possible  to  define  a  width  for  the  front  that  matches  the  observations.  For  the  scalar  wave 
equation  (Hedstrom  [1975]),  dimensional  analysis  reveals  the  width  of  the  front.  We  can 
not  use  the  same  technique  here.  Thus,  we  must  motivate  a  “definition”  of  width  of  front 
which  matches  what  we  observe  when  we  look  at  a  graph  of  the  solution.  We  will  show 
that  provides  us  with  a  definition  for  the  width  of  the  front. 

We  will  start  by  presenting  some  computations  of  the  path  of  steepest  descent  and  the 
height  of  that  path  for  a  particular  value  of  h.  These  computations  are  only  used  to  provide 
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motivation  for  both  a  definition  of  width  of  the  front  and  the  theorem  that  follows,  and  are 
not  used  for  proof.  However,  the  geometry  of  the  surface  is  complicated  enough  that 
a  specific  example  will  help  in  understanding  the  proof. 

Since  we  are  interested  in  the  behavior  in  the  vicinity  of  x  —  t,  we  look  at  contour 
plots  of  0?<£_  for  u  near  1.  From  the  topographical  interpretation  of  saddle  points,  we  define 
9?<£_  as  the  height  of  <p_,  since  \e*\  =  e®*. 


Figure  2:  Contour  plot  of  <£_  for  w  =  0.085  and  h  =  0.05.  Boldface 
lines  are  lines  of  steepest  descent  (constant  0<£_),  others  are  lines 
of  constant  height  (9?0_). 
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Figure  3:  Contour  plot  of  <t>_  for  w  =  0.D86  and  h  =  0.05.  Lines 
have  the  same  meaning  as  the  previous  figure. 


Figure  4:  Surface  plot  corresponding  to  figure  2.  The  direction 
of  view  is  from  the  upper  right  corner  of  figure  2  to  the  lower  left, 
viewed  from  above.  Note  the  branch  cut  in  the  center  of  the  plot. 

Note  how  the  path  of  steepest  descent  changes  between  figures  2  and  3.  The  path 
in  figure  2  has  three  components.  Two  parts  of  the  path  pass  through  the  saddle  points 
on  the  real  axis,  and  the  other  part  loops  around  the  branch  cut  at  j).  In  figure  3, 
the  path  of  steepest  descent  passes  only  through  the  principal  saddle  on  the  imaginary 
axis  (the  paths  through  the  other  saddle  points  are  forbidden  because  they  pass  on  the 
wrong  side  of  the  branch  cut).  What  we  are  seeing  here  is  a  complicated  interaction  of  four 
saddle  points.  Well  behind  the  front,  the  path  of  steepest  descent  has  two  features.  First  is 
the  part  of  the  path  that  mimics  the  path  for  the  continuous  problem.  This  is  the  part  that 
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loops  around  the  branch  cut.  The  other  part  consists  of  the  paths  that  pass  through  the 
two  saddle  points  on  the  real  line.  These  explain  the  oscillations  that  are  known  to  extend 
well  behind  the  front  (see  the  subsection,  ‘Behavior  Away  from  the  Front'). 

As  we  near  the  front,  die  paths  come  closer  together,  eventually  crossing  at  a  value  of 
a;  of  about  0.985.  It  is  interesting  to  look  at  the  height  of  the  path  in  this  critical  region.  If 
the  path  has  humps  at  each  saddle  point,  then  the  saddle  points  may  be  treated  separately. 
If  not,  then  we  must  consider  ail  the  saddle  points  together.  The  next  two  figures  answer 
this  question. 


i 
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and  the  height  of  the  path  for  w  =  0.985  and  h  =  0.05.  Note 
how  the  height  of  the  path  flattens  out  in  the  neighborhood  of  the 
saddles  ~6+,  -0_,  and  0_.  The  abscissa  in  the  second  graph  is  an 
arclength  along  the  path. 
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Figure  6:  Plot  of  the  left  side  of  the  path  from  the  principal  saddle 
and  the  height  of  the  path  for  u  =  0.986  and  h  =  0.05.  Note  how 
quickly  the  height  falls  as  the  path  moves  away  from  the  saddle 
points. 


These  figures  show  that  we  must  consider  all  of  the  saddle  points  simultaneously  in 
the  neighborhood  of  the  front.  However,  note  that  the  principal  saddle  point  6_  is  always 
higher  than  the  other  saddle  points.  This  means  that  the  principal  saddle  controls  the  gross 
features  of  the  solution.  In  particular,  we  expect  the  width  of  the  front  to  be  reflected  in 
the  behavior  of  the  principal  saddle.  The  next  figure  shows  the  behavior  of  the  height  of 
!  the  principal  saddle. 


Figure  7:  Height  of  the  principal  saddle  0_  as  a  function  of  w; 
h  =  0.05.  Note  how  rapidly  the  height  falls  as  w  increases  past  1. 

Since  5 R<£_(0_)  falls  off  so  quickly,  we  can  expect  that  by  considering  how  9ty_(0_) 
behaves  ahead  of  the  front,  we  should  be  able  to  find  the  width  of  the  front.  This  leads  to 
the  following  definition: 

Definition  1 :  Given  c  >  0  and  some  7,  >  0,  IF  the  height  of  the  principle  (highest)  saddle 
point  tends  to  -00  for  w  =  1  +  ch 1 ,  for  7  <  7,,  and  the  height  of  the  principle  saddle  is 
bounded  below  for  w  =  1  +  ch1,  for  7  >  7,,  THEN  the  width  of  the  front  is 

This  definition  is  based  on  the  following  estimate  of  the  inverse  transform  integral. 
Under  the  conditions  in  the  definition,  we  can  write  the  inverse  transform  integral  as 
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At-*™  f  ^ 


where  i /f(h)  =  o(l),  SR /(/»)  >  0,  and  £0  is  the  principal  saddle  point.  Then,  as  h  -*  0,  the 
integral  vanishes  (discarding  contributions  from  poles  near  the  path;  these  give  the  —  ^  in 
the  solution  ahead  of  the  front).  We  do  not  need  to  consider  the  effects  of  the  other  saddle 
points  because  of  the  exponential  importance  of  <j>  at  the  principal  saddle  point.  Note  that 
this  definition  does  not  say  that  the  width  of  the  front  is  asymptotically  equal  to  h ,  only 
that  the  width  of  the  front  is  at  most  hlc . 

The  restriction  in  definition  1  to  c  >  0  is  based  on  the  before-mentioned  difficulty  in 
defining  a  width  of  the  front.  For  c  <  0,  we  have  u>  <  1  and  hence  we  are  looking  behind 
the  front.  But  behind  the  front  there  is  no  clear  way  to  measure  the  width  of  the  front  (see 
the  graphs  of  the  computed  solution  in  section  5). 

For  more  detail  of  the  behavior  of  the  solution  near  the  front  we  must  take  into 
consideration  the  interaction  of  all  four  saddle  points.  This  is  discussed  below  in  the 
subsection  on  the  canonical  form.  We  proceed  below  to  find  the  width  of  the  front,  proving 
the  conjectures  above  about  the  behavior  of  the  path  and  the  height  of  the  principal  saddle. 

Theorem  2:  The  width  of  the  front  in  the  solution  to  (3.1)  is  0{h 2/s),  where  h  is  the  spatial 
step  size. 

Proof: 

First,  we  need  to  know  something  about  the  path  of  steepest  descent  (henceforth 
path).  We  actually  will  need  only  qualitative  information  on  the  behavior  of  the  path.  For 
this,  we  look  at  <£"(£)  for  £  a  saddle  point.  Now, 


2ih2(l  —  2  sin2  0)  +  8t  sin4  0 
*M)= - 


Kt  - ') 


Since  4  sin2(0)/fi2  —  l  is  negative  for  0  =  0_  and  positive  for  0  =  0+,  we  can  easily  evaluate 
<i>"_  at  the  saddle  points.  At  the  upper  saddle,  0_,  it  is  negative.  At  the  lower  saddle,  -0_,  it 
is  positive,  and  at  the  left  and  right  saddles,  T0+,  it  is  Ti~i,  for  some  real  positive  7.  These 
facts  allow  us  to  determine  the  direction  of  the  path  of  steepest  descent  at  the  saddle 
points;  see,  for  example,  Bleistein  and  Handelsman  [1975],  ch.  7. 
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Figure  8:  Saddle  points  of  <j>..  Arrows  point  downhill.  The  branch 
cut  is  indicated  by  the  shaded  line. 


Figure  8  shows  the  behavior  of  !R<A-  in  the  neighborhood  of  its  saddle  points.  Recall 
that  the  domain  is  the  strip  [ — 7r //i,  tt/Zx]  x  (-oo,oo]  in  the  £  plane.  Clearly,  the  path  of 
steepest  descent  is  one  of  three  forms  in  the  next  figure. 
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Figure  9:  Possible  paths  of  steepest  descent  for  <f>_.  Arrows  show 
the  direction  of  integration. 

We  first  study  the  possible  paths  of  steepest  descent  illustrated  in  the  previous  two 
figures.  Determining  which  of  these  three  choices  is  correct  for  all  values  of  h  and  u>  is 
very  difficult,  and  we  will  not  do  it.  Fortunately,  by  purely  geometric  considerations,  we  can 
show  that  the  saddle  point  at  0_  dominates  all  of  the  possible  paths. 

Recall  that  the  definition  of  the  height  of  a  saddle  point  at  0  is  W<£(0).  Figure  9a  shows 
the  path  which  passes  only  through  0_;  this  is  the  correct  path  if  9?<£_(0_)  is  less  than  the 
heights  of  the  right  and  left  saddle  points.  Of  course,  figure  9a  may  be  the  right  path  even 
if  this  is  not  the  case,  but  we  will  see  in  a  moment  that  that  is  not  important  to  us  here. 
In  the  other  cases  (figures  9b  and  9c),  the  saddle  at  0_  must  be  higher  than  the  other 
saddle  points,  for  otherwise  the  path  in  figure  9a  would  be  a  lower,  and  hence  better,  path. 
Therefore,  because  the  heights  of  the  saddles  enter  as  e*(e\  the  saddle  at  0_  is  always  the 
most  important  term.  We  will  call  this  saddle  point  the  principal  saddle  point.  This  is  net  to 
say  that  the  other  saddle  points  are  unimportant;  we  will  see  that  the  left  and  right  saddle 
points  lead  to  the  oscillations  that  are  observed  in  the  solution  of  (3.1).  The  lower  saddle 
point  —  0_  is  also  important  to  the  detailed  behavior  of  the  solution. 
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We  will  have  to  deform  the  path  of  integration  to  pass  through  the  principal  saddle 
point  before  we  proceed.  In  doing  so,  we  will  pick  up  a  term  from  the  singularities  along 
the  path  of  integration.  In  deforming  the  original  path  of  integration  passing  along  the 
real  axis  above  the  branch  cut  to  the  path  of  steepest  descent,  we  move  the  path  off  of 
singularities  at  0  and  ±  sin  '(/i/2)/2.  Only  the  pole  at  zero  contributes  anything.  As  with 
the  continuous  problem,  we  must  add  in  times  the  residue  of  this  pole  to  the  integral.  It 
is  easy  to  see  that  this  is  just  --2  by  making  a  Taylor  expansion  of  the  integrand.  We  will 
ignore  this  constant  term  in  what  follows. 

We  now  find  the  behavior  of  0_  and  4>.(0 _)  for  u  —  1  +  ch1,  as  functions  of  7  and  as 
h  goes  to  zero.  To  discover  what  sort  of  perturbation  expansion  we  will  need,  we  consider 
an  expansion  for  ui  =  1.  We  have  sin20_  =  -\h,  or  h£_  =  0_  «  isjhj'l.  Thus  we  see 
that  as  h  -*  0,  is  large  and  0_  is  small.  We  thus  proceed  by  making  expansions  with  £ 
large  and  0  small.  We  will  use  the  notation  .  •  •;  /»“)  to  mean  0(/tm,n^°' ' 

Using  this  result  as  a  guide,  we  let  w  =  1  +  ch1 ,  where  c  and  7  are  positive.  In 
addition,  we  will  assume  that  7  <  1.  It  is  easy  to  show  that  <£-(0-)  tends  to  —  as  h  — ►  0 
for  7  >  1. 

We  start  with  (3.6): 


•  2/1  1  - 
sin  0_  =  — - — 


w2~  1 


1  + 


w2h2 


(1-w2)2 


—  (1  +  u>)c/i7  —  (1  +  w)c/t7^ 

(I 

iin2  0_  =  — ch7^ 


2(1  +  w)2c2 

To  simplify  the  notation,  we  will  replace  (1  +  w)c  with  c.  Then  we  have 


)]' 


,  1  w2ft2-27 

1  +  T - o - 


+  0{h* 


or 

sin0_  =  iV c/i7^  1  +  ^-—2 - 

Since  0-  is  small,  we  may  expand  the  inverse  sine  to  get 


0 


1  w2/i2-2n 
He2) 


i\f  r:i  h2~< 
fj 


+  0{hA-7l^-h5l/2). 


Now,  we  can  find  0_(0.  ).  First,  we  need  to  know 


4  sin2  0_ 
h2 


-  1  =  -Ach1~2^\  +  )  +  0(/»2~37)  -  1 

-  I  =  2tVc/rr^^l  +  ^—--^(1  +  0(/»1-'l7;/»2-7)). 


(3.7) 
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Thus,  for  4>  (0  ),  we  plug  these  in  to  get 


-  ~2iy/ehi  2(l  + 


8  c‘ 

1  w2/i2  27 ' 


Iui2/i2  Zl\  iU 

+  H  ,,  ) 


+  ttV 

+  0(/i:J~7'r/2;  /*'  1/2;A5',/2  1 


Hcli1)2/2  h  1 
6 


Collecting  terms,  this  is 


*-(*- 


-■  4-  ( 1  —  a; 


Vch ~i~  2[  1  + 


(  1  uj2lt2  2l\ 

K -V-) 


wc-’/H’T'/8"1 

’  6 


+  0(/»3'7',/2;/»1_7/2;/»57/2'1) 

Now,  we  replace  I  —  w  by  —ch1  /[  1  +  u>)  (since  we  redefined  r  above)  and  get 


*-(*-) —2 
__  1 
~  _2 


—  c*/2#*37/*  ‘f— 0(/»3  77^2;  /il  ‘  7/2;  /ir'7^2  M 

V  I  +  w  6  / 

-  (positive  quantity)/i37',2~I  +  0(/i3  7l/2;h'  7/2;/ir>1/'2  ') 


(3.8) 


From  this  we  can  see  that  for  7  >  2/3,  <£_(0_)  -»  -  g  as  /»  — ►  ().  For  7  <  2/3,  0.(0 -)  -*■ 
-00  as  /j  -*  0.  Thus,  for  7  <  2/3,  U  -»  —  |  as  h  — *  0  (-  comes  from  the  pole; 
the  contribution  from  the  path  tends  to  zero,  as  does  the  term).  For  7  >  2/3  the 
contribution  from  the  path  does  not  tend  to  zero.  | 


Note  that  we  have  not  shown  that  (/  is  different  from  —  \  for  7  >  2/3,  since  we  would 
need  a  lower  bound  on  the  integral.  Thus  we  can  only  say  that  the  width  of  the  front 
is  0{h2/2)  rather  than  the  stronger  statement  that  the  width  is  asymptotically  proportional 
to  h 2/2 .  Section  5  will  present  computations  that  demonstrate  that  the  width  is  indeed 
proportional  to  h2/ 3.  Also,  this  result  is  good  for  t  small  and  positive  as  well,  since  the 
behavior  is  so  strong. 

The  limited  form  of  this  result  is  caused  by  the  complicated  behavior  of  the  saddle 
points  and  is  not  an  inherent  limitation  of  the  method.  With  better  information  about  the 
location  of  the  saddle  points  and  the  path  of  steepest  descent,  we  coulo  describe  the 
behavior  of  the  front  in  more  detail.  In  particular,  for  any  given  h  and  u>,  we  could  work 
out  the  detailed  behavior.  Since  we  are  interested  in  general  results,  however,  we  will  not 
do  this. 

This  sort  of  analysis  may  seem  contrary  to  the  analysis  in  section  2,  where  great 
care  was  exercised  to  get  the  contributions  from  all  of  the  saddle  points  along  the  path  of 
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steepest  descent.  This  is  because  this  section  and  the  previous  one  have  somewhat  different 
objectives.  In  the  previous  section,  it  was  important  both  to  introduce  the  techniques  of 
uniform  asymptotics  and  to  show  the  solution  to  the  continuous  problem.  Here,  we  are 
interested  in  finding  the  width  of  the  front  as  a  function  of  h.  To  compute  the  solution 
near  the  front,  we  would  have  to  determine  the  correct  path  of  steepest  descent  and  use  a 
uniform  expansion  to  calculate  the  behavior  near  the  front.  Using  the  canonical  form,  we 
do  this  computation  in  section  5. 

The  Canonical  Form 

In  section  2  we  saw  that  the  integral  representing  the  solution  of  the  continuous 
problem  could  be  transformed  into  an  integral  with  a  simpler  form.  We  did  this  by  simplifying 
the  argument  of  the  exponential  in  the  integral  into  a  form  which  represented  the  saddle 
point  behavior.  This  form  is  the  canonical  form,  and  by  studying  these  integrals  we  can 
gain  a  better  understanding  of  the  exact  solution  to  the  difference  equations. 

For  the  difference  approximation  of  this  section,  we  will  use  the  canonical  form  to 
discuss  the  behavior  of  the  principal  terms  in  the  asymptotic  expansion  near  the  front. 

The  canonical  form  chosen  depends  on  the  region  of  interest.  Near  the  front,  we 
know  from  above  that  9ty_(0_)  — j.  It  is  easy  to  see  that  9ty_(0+)  =  3ty_(-0+)  —  -5 

everywhere.  Therefore,  we  expect  the  path  of  steepest  descent  to  be  more  like  figure  9b 
or  9c  than  9a.  In  this  case,  we  will  need  to  consider  the  interaction  of  four  saddle  points. 
As  in  the  continuous  problem,  we  will  be  guided  to  a  canonical  form  by  considering  an 
expansion  of  <£_.  We  recall  that  the  saddle  points  are  Oj.  —  0(h l/2)  near  the  front.  Thus, 
in  contrast  to  the  derivation  in  section  two,  we  will  remain  in  the  0  plane  (rather  than  going 
to  the  l/£  plane). 

The  small  /»,  small  0  expansion  of  <t>_  (in  terms  of  ()  is 

•mo  ~  -  \ +  it + i(w  - l)* +  ^ (3-9) 

(In  this  expansion,  it  is  important  to  expand  the  sine  inside  the  square  root,  rather  than 
expanding  the  square  root  and  then  the  sin,  because  the  first  term  in  sin 20/h2  is  large, 
while  the  succeeding  terms  are  small.) 

Note  that  since  for  w  «  1,  (  se  and  hence  the  fc2£3  term  is  not  smaller  than 

the  first  two  terms.  However,  the  next  term,  £5h*,  is  smaller.  This  suggests  the  canonical 
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form 


ip(p)  —  ~  +  b  +  cp  +  dp3 4  =  <A-(0-  (3.10) 

P 

Here  we  have  already  applied  the  requirement  of  symmetry  of  the  saddle  points  to  eliminate 
the  p2  term.  By  setting  6  =  we  do  not  affect  the  saddle  points  of  ip.  We  may  also 
normalize  the  mapping  by  setting  one  of  the  remaining  parameters;  we  set  a  —  i/8  so 
that  d£/dp  se  1.  We  are  then  left  with  the  two  parameters  c  and  d.  To  determine  these, 
we  need  two  and  only  two  conditions.  These  two  conditions  are  provided  by  the  fact  that 
the  mapping  (3.10)  is  locally  one-to-one  at  the  saddle  points.  This  means  that  the  saddle 
points  correspond.  That  is,  if  f  is  a  saddle  point,  then  so  is  p(().  Since,  by  symmetry,  the 
saddle  points  are  described  by  two  parameters,  we  have  the  two  conditions  that  we  need 
to  determine  (3.10). 

Now,  ip(p)  has  four  saddle  points  p±  and  —  p±  which  satisfy 


-c  ± 


C2+  - 


3 id 
2~ 


6d 


(3.11) 


We  expect  (3.11)  to  have  two  real  and  two  imaginary  solutions.  This  will  be  the  case  if 
c2  +  3id/2  <  0.  There  are  many  ways  to  satisfy  this  condition;  it  is  clearly  satisfied  if  c 
and  d  are  pure  imaginary  and  d  is  small  enough  or  has  positive  imaginary  part.  This  is  just 
the  case  for  <p_.  Also  note  that  the  behavior  of  the  saddle  points  follows  that  of  </>_(£)  as 
h  -*  0.  In  particular,  the  saddle  points  of  ip[p)  go  toward  infinity  for  w  =  1  as  h  tends  to 
zero.  All  we  have  to  do  is  to  determine  c  and  d  such  that  ip(p±)  = 

Unfortunately,  it  is  too  difficult  to  analytically  determine  c  and  d.  We  can  solve  for 
them  numerically,  however,  and  compare  them  to  (3.9).  A  fixed  point  iteration  suggested 
by  Hedstrom  [1979]  for  the  scalar  wave  equation  can  be  adapted  for  this  problem: 

1  Pick  c(°*  and  d^°\  Set  n  =  0.  Define  V’^Hp)  by  (3.10),  with  and  d (")  replacing  c 
and  d. 

2  Repeat  steps  3-4  until  and  d^O  converge. 

3  Find  the  saddle  points  of  ip^.  Find  a  correspondence  between  these  saddle  points 
and  the  saddle  points  of  <p_. 

4  Solve  the  two  (linear)  equations  ip^n+l'l(p±)  =  <P-(Z±)  for  the  coefficients  c*n+1)  and 
d<n  +  ,>.  Set  n  +-  n+  l. 
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We  use  this  algorithm  to  find  c  and  d  for  w  in  a  neighborhood  of  l.  Good  choices  of 
c(0>  and  d{0)  are  needed.  We  can  start  at  w  =  1  by  taking  c  and  d  from  (3.9).  With  the 
u  =  1  values  of  c  and  d ,  we  can  proceed  to  nearby  values  of  w  by  continuation. 

This  canonical  form  defines  a  new  special  function,  which  is  a  generated  regular 
special  Bessel  function.  Graphs  of  this  function  are  displayed  in  section  5.  The  following 
two  figures  show  c/i  and  d/ih*  for  h  =  0.001.  The  figures  compare  these  to  the  values 
from  (3.9)  for  c  and  d. 


Figure  10:  Comparison  of  coefficients  for  model  equation  with  co¬ 
efficients  for  map.  Curved  line  is  c/i  for  h  =  0.001.  Straight  tine  is 
the  value  from  (3.9). 


d/(h**2  1) 


hop 


Figure  1 1:  Comparison  of  coefficients  for  model  equation  with  co¬ 
efficients  for  map.  Curved  line  is  d/ih2  for  h  =  0.001.  Horizontal 
line  (at  1/6)  is  the  value  from  (3.9). 


These  figures  show  that  near  the  front  the  form  in  (3.9)  (really  a  model  equation)  is  a 
good  approximation.  Away  from  the  front  the  neglected  terms  in  the  Taylor  series  in  (3.9) 
become  important,  and  (3.9)  is  no  longer  as  good  an  approximation.  Fortunately,  we  are 
interested  only  in  the  near-front  behavior,  so  we  will  use  the  (3.9)  form  of  the  mapping. 

We  can  use  the  above  result  to  look  at  the  behavior  of  the  canonical  form  near 
the  discontinuity.  Because  the  form  in  (3.9)  is  so  close  to  the  canonical  form  near  the 
discontinuity,  we  will  look  only  at  the  form  in  (3.9). 

Written  with  (3.9),  the  form  of  the  integral  for  the  solution  to  (3.1)  is 


+  i(<jj  -  1)  + 


We  now  consider  various  scalings  of  £  which  reveal  the  importance  of  the  terms  in  the 
exponential.  If  we  define  h°i)  =  £,  the  integral  is 


s'x,>  {(«/^  + i(“  - + 
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If  (3.1)  were  the  scalar  wave  equation,  a  =  -1  would  be  the  natural  scaling  (see  Hedstrom 
[1975]).  Writing  the  above  with  a  =  -1,  we  have 


/-.9“p{(w +i("“ l)’’ + 


(3.12) 


The  integral  (3.12)  is  useful  when  we  may  neglect  the  »7i2/8j?  term,  for  then  (3.12)  may 
be  written  in  terms  of  the  well  understood  Airy  functions.  For  that  term  to  be  negligible, 
we  must  have  |/»2/nl  <  |(w  -  1)t;|  and  \h2/r) |  <£  ]r?3|  at  the  saddle  points.  We  may  rewrite 
these  conditions  as 

|w  -  1|  » 


2 


(3.13) 


\V‘ 

|»?|  »  Vh. 

Now,  near  the  front  we  know  that  £  sa  l /'/h,  so  r\  =  —  >/h  and  the  condition  \r/\  >  \/h 

is  not  satisfied.  Thus  (3.12)  is  not  useful  near  the  discontinuity.  As  we  move  away  from 
the  front,  »?=/»$  ^constant  (from  (3.6)),  and  both  of  the  conditions  (3  13)  are  satisfied. 

In  particular,  (3.12)  and  (3.13)  show  that  the  term  in  the  integral  due  to  the  coupling 
(the  i7i2/8r;  term)  becomes  less  and  less  important  away  from  the  front.  From  this,  we  can 
conclude  that  in  a  region  near  (near  for  expansion  (3.9)  to  be  valid)  but  bounded  away  from 
the  discontinuity,  the  behavior  of  the  (/  term  is  similiar  to  the  behavior  of  the  centered 
difference  approximation  to  the  scalar  wave  equation. 

It  is  important  to  note  that  the  (/+  term  is  important  away  from  the  front  (Majda  an4 
Osher  [1977]),  so  we  may  not  exclude  it  in  discussing  the  behavior  of  the  solution  to  (31). 

Another  useful  scaling  of  £  is  a  =  —1/2.  With  this  scaling,  the  positions  of  the  saddle 
points  stay  nearly  fixed  in  the  7/  plane  (as  functions  of  h)  for  ui  near  1.  The  integral  then 
looks  like 

For  (3.14),  we  want  the  term  i(w  -  to  be  unimportant  and  t\fh  large.  For  this, 

we  have  |l/?/|  >  |(w  -  and  jr/3|  »  |(w  -  l)r///i|,  or 

|(w-l)| 


1 

In2! » 


h 

(w-  1) 


(3.15) 


From  (3.7)  and  rj  =  i/'/h ,  we  have  r?  «  hS1  '^2  at  w  —  1  =  ch1.  Inserting  these  into 
(3.15)  shows  (3.15)  satisfied  for  7  <  I  ((3.7)  is  valid  only  for  7  <  1).  Note  that  the  argument 
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of  the  exponential  (dropping  the  »(u  -  I  )/h  term)  has  4  saddle  points  equally  distributed 
about  the  unit  circle.  This  shows  us  that  very  near  the  discontinuity,  ail  four  saddle  points 
are  important. 

Note  also  that  this  scaling  makes  asymptotic  analysis  difficult,  as  we  would  need  iv/f> 
to  be  large  as  h  tends  to  zero. 


Behavior  away  from  the  Front 

We  may  use  the  integral  representation  (3  3)  to  9et  some  information  about  the  behavior 
of  the  solution  away  from  the  front  as  well.  Let  us  consider  the  effect  of  the  two  real  saddle 
points  in  figure  7.  There  is  some  *jc  <  I  such  that  for  w  <  ujci  the  two  saddle  points  on 
the  real  line  are  isolated',  the  path  of  steepest  descent  through  each  of  these  two  saddles 
starts  at  infinity  and  passes  through  only  one  of  them  before  going  back  out  to  infinity  (cf. 
figure  8c).  Thus,  we  may  use  classical  saddle  point  theory  to  estimate  the  contributions  of 
these  terms.  Considering  just  the  contributions  from  these  terms,  we  have 


Uoac  =  g(0+)~  +  g(-0+)  -  1 - +  ••• 

where  g  is  the  unexponentiated  part  of  the  integral  in  (3.3).  We  may  use  the  fact  that  0,  is 
real  to  relate  g(0+)  to  g(-0+)  and  <t>~{0+)  to  V  The  leading  term  in  g(0 ,)  is  h/0,  so 

to  leading  order,  g  is  odd  in  0.  We  also  have 

.  .  .  .  1  M  sin2  0+ 

*-(*»♦)  =  T  y  — >  ± 

Combining  these,  we  may  write  Uoae  as 


Uo,C 


VWio.)\ 


t/2 


sin 


4  sin2  0 
h2~ 


This  shows  where  the  oscillations  in  the  solution  to  (3.1)  come  from.  They  are  artifacts 
caused  by  the  presence  of  the  spurious  saddle  points  on  the  real  line  (spurious  in  the 
sense  that  they  are  not  present  in  the  continuous  problem  in  section  2). 

While  the  above  analysis  shows  where  the  oscillations  come  from,  it  does  not  give  the 
frequency  of  those  oscillations.  Since  0+  is  a  function  of  w  and  hence  of  x  and  t,  both 
g(0+)  and  ^'{0+)  contribute  to  the  oscillations.  Still,  this  does  show  where  the  oscillations 
in  the  solution  come  from,  and  a  more  careful  analysis,  such  as  those  in  Hedstrom  [1975] 
and  Pearson  [1969]  for  the  scalar  wave  equation,  would  show  the  observed  frequency. 
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Size  of  the  U+  term 

In  closing  this  section,  we  show  that  the  U+  term  in  (3.3)  is  asymptotically  negligible 
as  we  have  claimed.  For  the  region  away  from  the  iront  (w  bounded  away  from  1),  we  can 
use  the  same  basic  techniques  as  before.  First,  we  determine  the  path  of  steepest  descent. 
Second,  we  determine  the  appropriate  canonical  form.  The  canonical  form  near  the  front 
will  be  different  from  the  canonical  form  away  from  the  front.  Finally,  we  get  estimates  on 
the  size  of  the  terms  in  the  resulting  integrals.  For  the  region  near  the  front,  the  situation 
is  more  complicated.  In  this  region,  we  will  show  that  the  U+  term  is  of  the  order  of  the 
first  neglected  term  in  the  t/_  expansion. 

First,  for  the  region  away  from  the  front,  we  look  at  the  saddle  points.  iscover 
the  path  of  steepest  descent,  we  consider  the  relation  between  <t>-{£±)  and  </>+(£*)•  Let 
0  =  0  +  ic.  Clearly, 


=<t>-(0)  4- 

h 


since  the  sign  of  the  sine  changes  the  sign  of  the  square  root.  Hence,  5Rc£+(0)  = 

Thus  the  behavior  of  the  saddle  points  and  the  path  of  steepest  descent  is  the  same  as 
that  for  4>_  found  earlier  in  this  section. 

There  are  then  two  cases  to  consider:  w  away  from  1,  where  the  saddle  points  are 
well  separated,  and  w  near  1,  where  the  saddle  points  are  coming  together  (at  infinity).  We 
will  first  consider  the  simpler  case  of  w  away  from  1 . 

In  this  region  of  u,  we  can  use  classical  saddle  point  analysis  to  get  an  estimate  on 
U+.  Following  Erdelyi  [1956]  we  have  that 

U+  ss  — —  *  e^+(*-)«(7  +  asymptotically  smaller  terms. 

I 

Here,  g(0_)  —  a(0_)(l  -f  A+)  (see  (3.3)).  We  can  investigate  the  small  h  behavior  of  this  by 
making  Taylor  expansions  in  h  of  g  and  <j>".  We  already  know  that  3? <f>+(9-)  «  —  5.  The 
first  terms  in  the  expansions  for  g  and  4>"  are  ()»{  =  0,  0  =  0  +  7r,  £  is  0(1)  away  from  the 
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front) 


9(0) 

K(0) 


/»2c(t  +  2 it  -  »v^2  -  r) 

8v/4£2  -  1 

2t 

(4£2  _  i)3/2 


Thus,  U+  is  0(h 2).  This  is  just  what  we  would  expect,  since  the  method  that  we  are 
considering  is  second  order  accurate.  Note  also  that  the  exponential  is  like 


=slowly  varying  part  X  e'x"/h 

so  we  would  expect  IJ+  to  have  a  wavelength  of  2h  in  x.  These  oscillations  are  observed 
in  practice. 

As  a  demonstration  of  this  analysis,  the  following  table  presents  the  results  of  computing 
U f  numerically  for  h  —  2-n,  n  =  5, 6, . . . ,  1 5. 


n 

u+ 

ratio 

5 

1.44E-2 

NA 

6 

3.88E-3 

1 .89727 

7 

9.49E-4 

2.03098 

8 

2.32E-4 

2.03180 

9 

5.71  E-5 

2.02369 

10 

1.41E-5 

2.01370 

11 

3.52E-6 

2.00731 

12 

8.77E-7 

2.00377 

13 

2.19E-7 

2.00191 

14 

5.47E-8 

2.00096 

15 

1.37E-8 

2.00048 

Table  1:  U+  for  h  =  ■2~n  and  x  =  2 1  —  5.  Ratio  column  is  the 
log  of  the  ratios  of  success  e  rows  divided  by  log(1/2).  Note  how 
quickly  U+  approaches  0(h 2).  The  value  for  U+  was  taken  as  the 
maximum  over  two  periods  (each  period  is  2 h  long)  and  a  constant 
factor  was  divided  out.  Of  course,  there  is  no  ratio  for  the  first 
entry. 

It  is  instructive  to  note  that  as  w  goes  to  1,  <t>"  becomes  0(/i3/2).  Thus,  as  w  nears  1, 
it  is  no  longer  sufficient  to  consider  just  the  effect  of  the  principal  saddle  point.  We  look 
then  at  the  canonical  form  (3.10).  By  the  above  analysis,  <j>+(0_)  and  <t>-(0_)  are  essentially 


the  same  (they  differ  by  a  term  constant  in  0).  Thus  we  can  use  the  same  canonical  form 
for  estimating  the  size  of  th  f/f  term. 

The  behavior  of  the  canonical  form  integral  is  very  complicated  for  w  «  1.  The  best 
that  we  can  do  with  it  is  to  show  that  it  is  comparable  to  the  first  neglected  term  in  the 
previous  analysis.  Define  to  be  the  unexponentiated  part  of  the  f/+  integral  in  (3.3) 
and  </_  to  be  the  unexponentiated  part  of  the  U .  integral.  Because  of  the  behavior  of  the 
saddle  points  of  4>+  and  <t>  near  the  front,  we  need  to  know  the  behavior  of  yt(0)  for  0 
near  n  and  the  behavior  of  y  (0)  for  0  near  0.  Let  0  =  it  +  0,  0  small.  This  gives  us: 


Now,  f(0)  =  i7icot(0/2)/2,  so  /(0  +  rr)  =  —ih  tan(0/2)/2.  From  this  we  can  see  that 
gt(0)  is  smaller  than  3  (0),  since  both  are  /  X  similiar  functions,  and  tan  0  <£  cot0  for  0 
small. 

We  have  not  shown  that  this  term  is  small;  indeed  we  will  not.  It  is  a  recognized 
problem  in  asymptotic  analysis  that  estimating  the  higher  order  terms  in  an  expansion 
for  arbitrary  h  and  u  can  be  very  difficult.  For  specific  h  and  w,  however,  these  terms 
could  be  estimated  by  preforming  the  integration  numerically  with  a  method  with  a  known 
error  bound.  See  Wong  [1980]  for  a  clear  discussion  of  error  estimates  for  asymptotic 
expansions  of  integrals. 

To  get  a  good  estimate  of  f/t ,  we  would  need  to  know  both  the  location  of  the  saddle 
points  and  the  path  of  steepest  descent.  Since  both  the  computations  of  the  path  of 
steepest  descent  and  the  integral  (3.14)  show  that  all  of  the  saddle  points  are  important, 
we  must  be  able  to  estimate  the  canonical  form  integral  corresponding  to  these  features. 
As  the  canonical  form  integral  represents  a  new  special  function,  this  analysis  reduces  to 
tabulating  the  integral  which  the  canonical  form  (3.10)  represents  and  using  the  resulting 
tables  to  estimate  f/+ . 

We  can,  however,  provide  computations  to  show  that  the  term  is  in  fact  as  inconse¬ 
quential  as  claimed.  For  this  specific  problem,  we  can  look  at  the  following  table: 
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mssmon 


n 

u+ 

ratio 

5 

2.02E-4 

NA 

6 

9.77E-5 

1.04968 

7 

2.57E-5 

1.92524 

8 

8.09E-6 

1 .66781 

g 

4.34E-6 

0.89907 

10 

1.61E-6 

1 .42923 

ii 

5.09E-7 

1.66131 

12 

2.26E-7 

1.17489 

13 

8.16E-8 

1.46689 

14 

3.19E-8 

1.35582 

15 

1.28E-8 

1.31484 

Table  2:  U+  for  h  =  2-n,  <  =  1 ,  and  w  =  1  +  A2/3.  Columns 
have  the  same  meaning  as  in  the  previous  table.  The  lack  of  an 
obvious  limit  in  the  ratio  column  is  due  to  the  need  to  average  over 
the  oscillations.  Unfortunately,  the  predicted  ratio  changes  over  the 
length  of  one  oscillation,  thereby  contaminating  the  results. 

For  more  general  problems  (such  as  the  next  section),  we  can  look  at  the  graphs  of 
the  appropriate  canonical  form  which  are  presented  in  section  5. 
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4.  Behavior  of  Difference  Approximations  near  the  Discontinuity: 

General  Difference  Schemes 

In  this  section  we  analyze  a  model  equation  for  the  model  problem.  This  approach 
was  discussed  in  section  1.  Our  model  equation  is 


du 

dv 

dp+'v 

,,adi+ lv 

+  ah? - - 

+  bh 9 - - 

dt 

dx 

dxp+l 

dxi+x 

dv 

du 

,„dp+lu 

+  ahp - r 

+  bhq - : — 

dt 

dx 

dxP^1 

dxi+l 

t 

4’ 

if  x  <  0 

:.0) 

=  v(a:,0) 

=  <  0, 

if  x  =  0 

1 

4  * 

if  x  >  0 

l  <P  <  q. 

Here,  p  and  q  are  integers,  and  one  of  p  and  q  is  even,  the  other  is  odd.  For  example,  for 
p  —  2,  a  =  -1/6,  p  —  1,  and  6  =  0,  this  equation  models  the  nondissipative  second  order 
scheme  studied  in  section  3.  We  will  use  p  to  allow  us  to  see  what  effect  the  lower  order 
term  has  near  the  discontinuity.  For  p  =  0,  (4.1)  is  just  the  model  equation  for  difference 
approximations  to  the  ordinary  wave  equation;  p  =  1  is  the  telegrapher’s  equation,  and  it 
is  this  case  that  should  be  kept  in  mind  for  most  of  the  analysis. 

That  this  model  equation  is  a  good  model  can  be  proved  by  techniques  similar  to  those 
in  section  3.  In  fact,  section  3  can  be  taken  as  proof  of  this  model  equation  for  the  case 
of  the  second  order  difference  approximation  in  (3-1). 

Averaging  of  the  undifferentiated  term  has  not  been  included  in  order  to  simplify  the 
analysis.  The  techniques  used  here  can  be  used  to  find  the  behavior  of  difference  schemes 
which  average  the  undifferentiated  terms. 

Our  analysis  will  proceed  along  the  same  lines  as  in  the  previous  section.  We  represent 
the  solution  to  (4.1)  as  an  inverse  Fourier  integral.  We  then  study  the  behavior  of  the  saddle 
points  of  this  integral,  and  show  that  the  width  of  the  front  is  The  calculations 

for  this  are  very  similiar  to  those  in  section  3  and  we  will  compare  the  results  there  to 
the  results  in  this  section.  We  then  use  p  to  show  in  which  regions  near  the  front  the 
undifferentiated  term  is  important.  Finally,  we  present  a  generalization  of  (4.1)  which  shows 
exactly  what  dimensionless  quantity  must  be  small  for  our  analysis. 


Integral  Form  of  the  Solution 


We  start  with  the  Fourier  transform  of  (4.1): 


du 

dt 

dv 

dt 


=  —  i£v  +  ai£(ih£)pi'  +  bi£(ili£)‘lv 
=  —  i£u  +  ai£(ih£)pii  I  bi£(ih£)vu  —  pv 


where  £  is  the  dual  variable.  In  matrix  form,  this  is 

±(U\  =  (  o  -*e( i  -  «(*/*€)*>  - 

dt\v)  V-*^1  -  a(ih£Y  -  b(ih£)q)  -p 

Then  the  eigenvalues  of  this  matrix  are 


/'(MOV 


(4.2) 


\±  =  -~±  j\Z't(2Yl  -  a(iA()P-  b(ih£)«)2  -  p2 

Let  <t>k  =  X±  +  iu£,  where  u  =  x/t.  The  branch  cuts  (there  are  many  of  them)  are  chosen 
so  that  the  h  =  0  limit  has  saddle  points  with  the  same  behavior  as  for  the  continuous 
problem.  The  other  branch  cuts  are  chosen  to  be  well  away  from  the  origin  (the  region  of 
interest)  for  u  sa  1.  If  we  let  a±  be  the  unnormalized  eigenvectors  of  the  matrix,  we  can 
write  the  solution  to  (4.2)  as 


(*)  =  «r(€K«x,‘  +  «-(0 


where  n±(£)  are  determined  by  the  Fourier  transform  of  the  initial  data.  Then  the  solution 
of  (4.1)  is  the  inverse  transform 


1  f+°° 

u(x,t)  =  —  y  [r*,(0«u',X  +  t  +  «  .(£)a 


i  _«x  - 


(43) 


where  «t±  are  the  first  components  of  a±. 

It  is  possible  to  express  «±(£)  in  terms  of  the  initial  data  and  Xt,  but  it  shows  nothing 
new.  We  will  therefore  assume  that  «_(()  \/‘l£  for  £  near  0,  and  f".<.  the  X+  term  is 

negligible  compared  to  the  X_  term  for  x.  «  t.  The  arguments  for  these  assertions  follow 
the  same  pattern  as  those  in  the  previous  section,  and  they  will  not  be  repeated  here. 


Width  of  the  Front 

To  determine  the  asymptotic  behavior  of  the  width  of  the  front,  we  will  follow  the 
procedure  used  in  the  previous  section.  Specifically,  we  first  find  the  saddle  points  of 
4>4£).  We  then  show  that  the  saddle  point  corresponding  to  the  principal  saddle  of  the 


continuous  problem  controls  the  width  of  the  front  Finally,  we  show  flow  <fi  (^(to  +  ch1)) 
behaves  as  a  function  of  7. 

Consistency  of  the  difference  approximation  will  be  used  to  guide  and  justify  the  steps 
in  the  analysis.  We  require  that  the  h  >  0  limit  of  the  equations  below  must  give  the 
respective  equations  from  the  continuous  problem.  In  particular,  we  will  use  our  knowledge 
of  the  behavior  of  the  principal  saddle  point  to  pick  an  appropriate  perturbation  expansion. 

Theorem  3:  The  width  of  the  front,  as  defined  by  definition  3.1,  under  a  general  difference 
scheme  with  model  equation  (4.1)  is  fip/(pt 

Proof: 

The  saddle  points  satisfy  the  equation 

„  ,  .„*(>  -  -  KihW) (I  -  (I  +  p)a(ih£)p  -  (1  +  <,)b(ihm  ,  , 

U  —  ±ll - —  _ .:r_ zz= - h  lu).  (4.4) 

\/'U2(l  -  a(i/iOp  -  b(ih()<> j2  -  p2 

Now,  from  section  2  we  expect  that  while  (  tends  to  infinity  as  w  approaches  1,  /*£  will 
remain  bounded.  We  define  a  new  variable  7/  =  ih(,  including  the  i  to  simplify  the  notation. 
We  are  now  looking  for  solutions  q  to  (4.4),  where  ?/  is  small  for  w  near  1.  Expanding 
(4.4)  as  a  polynomial  by  moving  the  iui  term  to  the  left  hand  side  and  then  squaring  both 
sides,  we  can  now  drop  the  terms  in  the  polynomial  from  (4.4)  with  r/  to  powers  greater 
than  < 1  +  2.  This  leaves  us  with: 

8(2  +  p  -  w2)«t?p+2  +  8(2  +  q-  w  2)6r;‘,+z  +  1r;2(w2  -  1 )  +  /r2w2p2  =  o(t7<7+2).  (4.5) 

This  equation  is  clearly  too  complicated  to  solve  analytically.  We  thus  turn  to  perturbation 
techniques.  Our  interest  is  in  h  «  0  and  1.  Starting  with  a  perturbation  series  in  h 
for  the  solutions  t)  to  (4.5),  we  have  as  the  equation  for  the  zeroth  order  term 

8(2  +  p  -  w2)«r?o  f  p  +  8(2  +  p  -  aj2)br]Q+q  +  i?7o(w2  —  1 )  =  0. 


The  first  p  +  2  roots  of  this  are 
rO  twice 

770  I  p  roots  of  t/ 


_2-cv2 
2a(2  +  p  -  w2  j 


t. 


I  —  c 0 

«( 1  +  p) 


+  higher  order  terms 


The  remaining  q  -  p  roots  we  may  ignore  because  they  do  not  approach  zero  as  w  -*  I . 
We  can  make  an  identification  of  these  saddle  points  by  comparing  with  the  scalar 


wave  equation  case  (Serdjukova  [1971]  and  Hedstrom  [1975]).  In  that  problem,  there  are 
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p  saddle  points  due  to  the  difference  approximation  about  the  origin;  they  correspond,  to 
first  order,  to  the  p  saddle  points  found  above.  The  other  two  saddle  points  correspond  to 
the  two  saddle  points  in  the  continuous  problem  discussed  in  section  two.  We  will  need  to 
know  more  precisely  the  location  of  these  last  two  saddle  points,  as  we  use  the  behavior 
of  the  height  of  the  principal  saddle  point  as  our  definition  of  width  of  front.  To  get  the 
next  term  in  the  perturbation  expansion  for  q  q0  =  0,  we  let  r/  --  -dha .  Here,  d  has 
been  chosen  real  and  positive  (so  that  £  =  idhai  will  be  the  principal  saddle  point).  Also, 
since  we  are  interested  in  the  width  of  the  front,  we  need  to  know  how  r\  changes  as  a 
function  of  cj,  with  u  =  1  +  ch1 .  Making  both  of  these  substitutions  in  (4.4),  we  get 

8(2  +p-  w2)a(-d)2+p/i(2+'pl“  4-  8(2  4-  q  -  .u2)fe(-d)<l42h(2  +  ,,)a 

+  4(w  4-  l)r(-d)2/i2Q  +  "’  +  /i2u;V  =  0. 

Now,  since  q  >  p,  we  may  drop  the  term.  This  leaves  us  with 

8(2  4-  p-  w2)a(-d)2+p/i(2+p!°  +  4(u>  4-  l)c(— d)2h2a4 11  +  h2tj2p2  =  0.  (4.6) 

There  are  three  cases. 


case 

condition 

exponents 

1 

2a  -f  7  <  2 

(2  4-  p)a  =  2a  4-  7 

2 

2a  +  7  >  2 

(2  4-  p)a  =  2 

2 

3 

2a  +  7  =  2 

all  equal 

2 

a  J  —d 

p/  (u»  4 

V 

VI 2  +  P)  I 


iils 


_  'izi 

2 


Table  3:  Three  cases  for  q.  The  condition  column  gives  the 
condition  for  that  case  to  hold.  The  terms  column  gives  the  powers 
of  h  in  (4.6)  for  that  case,  and  the  last  two  columns  give  n  and  — d 
for  that  case.  The  *  value  of  — d  is  the  solution  to  a  +  2-tn  order 
polynomial,  and  is  omitted. 

Each  of  the  cases  in  table  3  is  valid  for  a  range  on  7,  or,  equivalently,  for  a  region  near 
the  front.  Case  1  is  valid  for  7  <  2/(1  4-  2/p).  Since  this  bounds  7  from  above,  the  region 
where  case  1  is  valid  is  asymptotically  further  from  the  front  than  cases  2  or  3.  Note  from 
(4.6)  that  in  case  1 ,  the  fi2w2p2  term  may  be  dropped,  and  hence  this  is  the  ■•>gion  where 
the  undifferentiated  part  of  (4.1)  does  not  affect  the  solution  (to  first  order)  Case  2  is  valid 


for  7  >  2/(1  +  2/p).  Since  this  bounds  7  from  below,  this  region  is  asymptotically  closer 
to  the  front  than  the  region  in  case  1.  In  this  case,  the  4(iu  +  l)e(—  ^7)2/l2Q+',  term  in  (4.6) 
(the  term  from  the  convective  or  z-derivative)  is  less  important  than  the  p  term.  We  note 
that  this  case  is  for  w  very  near  the  discontinuity,  closer  than  we  expect  the  approximation 
to  the  front  to  be.  Case  3  is  the  transition  region  between  cases  1  and  2.  In  this  region, 
all  parts  of  (4.1)  are  equally  important. 

Now  that  we  know  something  about  the  behavior  of  the  principal  saddle  points  for 
u  ra  l,  we  can  look  at  the  path  of  steepest  descent.  As  in  the  previous  section,  we  will 
need  very  little  information  about  the  path  in  order  to  find  the  width  of  the  front.  The 
following  figure  shows  the  situation  for  the  methods  in  this  section. 


Figure  12:  Possible  paths  of  steepest  descent  for  <f>_  near  the 
saddle  points.  Arrows  show  the  direction  of  integration. 

As  before,  the  path  always  passes  through  the  principal  saddle  point. 

We  now  look  at  the  value  of  <f>_  at  the  principal  saddle  point  r;  =  -dh°.  It  will  be 
useful  to  expand  <p_  at  this  point  for  £  large  and  r/  small 


50 


<t>  (0  —  -  *  -  *  n/^0  -  *  *  *<v£ 

=  -  *  -i^O  —  />( t/<£),»)2  ^  4-  j'w^. 

Expanding  the  argument  of  the  square  root  we  get 

<M0  =  -  2  -  1  -  -(l,lp  -  2bri'>  -  --*2  4-  C(r/2P)  4-  iw£ 

—  ~  ~  a,'P  ~  b'lq  “  + 

Collecting  terms,  we  have  finally 

4>  (0  =  ~  *£0  ~  w)  4-  zf^ar/P  4-  tr?’  4-  ^  4-  0(72p)j 

or 

<t>  (0  —  *£0  ~  u  -  arip)  +  r/2p) 

Now  we  substitute  —i?i/h  for  £  and  —ch i7'for  I  —  w  to  get 

<t>  (0  =  -  *//»'  1  (—ch.1  -  a//p)  +  0(r),+  l/t  1 ;  /i2tj-2; 7]2‘l). 

Finally,  we  have  y  —  -dh°,  or 

<M0  =  —  dka  '[\ch1  +  a(—d)p  hJ><*)  4-  0(/ta(,+  ,)_1;  h2~2°)  h2pa). 

Since  d  >  0,  we  clearly  need  only  look  at  the  behavior  of  ha  l(ch’1  4-  c(-d)phpa).  We 
wish  to  find  the  critical  7C  such  that  this  term  (which  controls  the  height  of  the  principal 
saddle  point  for  w  »  I)  approaches  -oo  for  7  <  7C  and  0  for  7  >  7 c.  Flowever,  we  can 
get  no  further  without  assuming  something  about  7.  We  will  first  take  case  1  from  above, 
in  this  case,  pa  —  7  so  both  terms  are  equally  important.  We  then  have 

■M0  =  -^—55)  +  -  («> 

where  we  have  used  the  formula  for  —  d  above.  The  quantity  in  parenthesis  is  always 
positive,  so  this  term  has  the  critical  gamma  7,.  =  l  -  a.  From  table  3,  we  see  that  case  1 
implies  that  (2  4-  7)^  2r*  4-  7,  so  we  have  finally  7e  =  p/(p  4-  I). 

We  need  only  show  that  the  other  two  cases  in  table  3  are  excluded.  First,  we  have 


pa  >  7  in  cases  2  and  3.  Since  we  are  looking  for  the  critical  7,  this  means  that  we 
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must  have  7  =  1  —  a.  From  table  3  we  know  that  in  cases  2  and  3,  a  =  2/(2  +  p),  so 
7  =  1/(2  +  p)  and  hence  that  2a  +  7  <  2.  But  for  both  cases  2  and  3,  the  condition  on  7 
is  2a  -t-  7  >  2.  Hence  for  7  =  -ye,  case  t  is  the  only  case  that  applies.  | 

In  the  proof,  we  again  expanded  about  h(  =  0.  (Recall  that  in  the  continuous  problem 
we  expanded  about  1/f  =  0.)  If  we  could  not  have  done  this,  then  we  would  not  have 
been  able  to  use  the  model  equation  (4.1),  since  it  is  essentially  an  expanssion  in  h£  for  a 
difference  approximation  to  (2.1). 


Results  of  Section  3  as  an  Example 

To  illustrate  the  results  above,  we  will  briefly  sketch  out  the  analysis  for  the  difference 
method  analyzed  in  section  3.  The  parameters  for  the  centered  second  order  difference 
method  considered  there  are  p  =  2,  a  =  -1/6,  p  =  1,  and  6  =  0.  Then  for  case  1  in 
table  3,  we  have 


-d  — 


-c(u  +  1) 


2(  —  1  /6)(2  +  2  —  w2) 


\/(l  +  w)e 


since  w2  ca  1 .  With  this,  we  may  write  (4.7)  as 


*-(«)  -  -5  -  w + -  s) + ' ' " 

(Recall  that  c(l  +  w)  here  is  c  in  section  3.)  This  is  the  same  as  (3.8).  Finally,  case  1  is 
valid  for  7  <  2/2  =  1,  and  the  width  of  the  front  is  0(/»2/3). 


The  Canonical  Form 

As  in  the  previous  section,  it  will  be  very  useful  to  construct  the  canonical  form  for  the 
solution  (4.3)  to  (4.1).  By  producing  graphs  of  the  canonical  form,  we  will  be  able  to  show 
the  near-front  behavior  for  a  wide  range  of  difference  approximations. 

Near  the  front,  we  have  the  following  small  h,  small  tj  expansion  for  <£_(£) 

<M0  «  ~  \  -  *€(1  -  w)  +  ~  +  ai£{ih£)p  +  biS(ihZ)'. 

Following  the  example  in  section  3,  we  expect  the  canonical  form  to  be 
ij){r)  —  -  -F  6  +  cr  -F  •  •  •  4-  drp+l  +  er,+  1. 

The  parameters  may  be  determined  by  the  method  outlined  in  section  3.  However,  the 
exact  form  of  ip  is  less  interesting  than  the  simpler  approximation 

ip{r)  =  -  i  -  t>(l  -  w)  +  +  air(ihr)p  +  bir(ihr )* 

2  8  r 


h 


* 


which  is  just  the  leading  terms  in  <£_(r).  It  is  straightforward  to  show  that  this  form  has  the 
required  properties;  p  +  2  saddle  points  which  converge  to  a  saddle  point  of  infinite  order 
as  h  — *•  0  and  u  — ►  1.  The  q  —  p  saddle  points  are  much  further  from  the  origin  than  the 
first  p  -i-  1  saddle  points  and  are  unimportant  to  first  order. 

Because  of  the  complexity  of  even  the  approximate  <f>{r),  we  will  not  analyze  it  here. 
Instead,  the  reader  is  refered  to  section  5,  where  graphs  of  the  canonical  form  integral 
are  presented.  These  graphs  show  both  the  behavior  of  difference  schemes  that  fit  the 
framework  of  (4.1)  near  the  front,  and  the  qualitative  similiarity  between  the  finite  difference 
solutions  to  the  scalar  wave  equation  ut  =  ux  and  the  telegrapher’s  equation. 


The  Cell  Wave^Coupling  Number 
Consider  the  generalization  to  (4.1) 


du 

dt 

dv 

dt 


4-  ahp 


dp^v 

dxP+' 


d u  ,  „ 

— c- — h  ahp 
dx 


dp+ 1  u 

di~p+i 


+  bhq 


di"v 
<9z9+  1 


,,„dq+lv 


(48) 


All  we  have  done  here  is  to  add  a  propagation  speed  c  to  (4.1).  The  advantage  of  this  is 
that  now  (4.8)  is  in  dimensional  form  if  a,  b,  and  c  have  the  dimensions  of  a  speed  and  p 
has  dimension  l/t.  Now,  we  chose  a  change  of  variable  to  put  (4.8)  into  nondimensional 
form.  Let  z  =  ex' /p  and  t  =  t' / p.  Define  It  —  hp/c.  We  will  call  this  the  cell  (from  the  h) 
wave-coupling  (from  the  p)  number.  Then  (4.8)  may  be  written  as 


du  dv  a,„dp+lv  b„adq+iv 

—  —  ___  +  - !{? - -t_  -Hq - 

dt'  dx'  c  Qx'P+i  c  dx'q+x 

dv  du  a  .„.dp+lu  b„adq+lu 

__  —  —  —  H — It  - -  +  -ltq~ - - —  v. 

dt'  dx'  c  dx,P+\  c  dx",+  i 

This  shows  that  the  appropriate  dimensionless  quantity  which  must  be  small  in  our  analysis 
is  It,  the  cell  wave  coupling  number.  This  is  illuminating,  as  it  shows  that  if  the  cellsize 
is  small  or  the  coupling  coefficient  is  small,  then  the  solution  to  (4.8)  will  behave  near  the 
front  like  the  solution  to  the  scalar  wave  equation  (p  =  0  in  (4.1)). 


5.  Numerical  Results 


We  present  in  this  section  two  distinct  numerical  results.  The  first  set  of  computations 
is  used  to  illustrate  the  theory.  They  show  that  the  width  of  the  front  is  as  predicted  for  a 
number  of  schemes,  both  semi  and  fully  discrete.  The  second  set  of  computations  presents 
graphs  of  the  generalized  Bessel  functions  described  in  section  4.  These  show  that  the 
qualitative  behavior  of  approximations  to  the  telegraph  equation  for  small  phjc  (recall  that 
this  is  the  cell  wave  coupling  number)  is  the  same  as  the  behavior  for  the  scalar  wave 
equation. 


Behavior  near  the  Discontinuity  for  Various  Difference  Schemes 

The  following  graphs  show  the  width  of  the  front  for  two  method-of-lines  (MOL) 
schemes  and  two  fully  discrete  schemes.  The  MOL  schemes  are  the  centered  second  (5.1) 
and  fourth  order  (5.2)  schemes 


dUi  _  Vj+j-Vi-i 
dt  2h 

dVj  =  Uj+t-Uj-t 

dt  2h 


(5.1) 


and 


dUi  =  +  i  -  Vj-i)  -  [Vi+2  ~  Vi- 2) 

dt  I2h 

dV±  =  _  —  Uj~i)  —  (Uj+2  —  Uj-2.)  _ 

dt  12 h  ' 


(5.2) 


for  t  =  0,  ±1,  ±2, . . .  and  h  the  space  step  size.  The  fully  discrete  schemes  are  Leapfrog 
(5.3)  and  Lax-Wendroff  (5.4).  Let  k  and  h  be  the  time  and  space  step  size  respectively, 
and  let  X  =  k/h.  Then  the  schemes  are 


and 


C/”+1 

yn+i 

w  t 


=ur'~  x(K?+.-vr_.) 


(5.3) 


u?+l  =4X(xi/7+,  +  uk  -  i)k?+1) + (i  -  x2K/r  +  +  (1  - 

v?+l  =^X(XK?+1+(i-A-  lK7,n+.)  +  (l  -*-X2+i*2)Vr  (5.4) 

+  *X(XK?_, +(!-**)(/?_,) 
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Each  method  was  integrated  to  t  =  1/4  for  h  s=  1/2"  for  n  =  5, . . .,  10.  For  each 
value  of  h,  the  width  was  determined  by  finding  the  value  a  x  where  the  solution  V  was 
equal  to  U0 .  The  values  of  U0  used  were  .05  and  .001.  These  were  chosen  to  measure  the 
width  in  two  different  areas:  inside  the  front  ( Uq  =  .05)  and  in  the  tail  ahead  of  the  front 
(Uo  =  .001).  By  chosing  two  different  values  of  U0,  we  show  that  the  width  behaves  as 
predicted  over  the  entire  front.  Since  U  is  defined  only  at  discrete  points,  we  used  linear 
interpolation  to  find  the  value  of  x. 

First,  we  show  the  computed  solution  for  h  =  1/512  and  t  =  1/4  for  each  method. 
The  MOL  schemes  were  integrated  with  RKF45  (see  Shampine,  Watts,  and  Davenport  [19/6] 
for  a  discussion  of  the  merits  of  RKF45;  Forsythe,  Malcolm,  and  Moler  [1977]  for  the  code); 
the  fully  discrete  schemes  used  k  =  h/ 1.25.  Following  graphs  of  the  computed  solution 
are  eight  graphs  showing  the  width  of  the  front  for  these  four  schemes. 
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Width  For  Lax-WendroFF 


u  =  .05 
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The  graphs  of  the  width  of  the  fronts  verify  the  predictions  of  the  theory  developed  in 
sections  3  and  4.  The  second  order  method  of  lines  calculations  verify  the  results  in  section 
3;  the  others  confirm  the  results  in  section  4.  Notice  that  for  the  fourth  order  method  of 
lines  the  width  of  the  front  is  as  expected.  Note  also  that  for  the  finite  difference 

methods  (Lax-Wendroff  and  Leapfrog)  the  rate  of  convergence  to  h 2/<3  is  slower  than  for 
the  semidiscrete  approximations.  This  points  out  the  asymptotic  nature  of  the  theory:  the 
results  are  applicable  only  for  large  t/h. 


Graphs  of  Generalized  Bessel  Functions 

We  have  seen  in  previous  sections  that  the  behavior  of  difference  approximations  to 
the  telegrapher's  equation  may  be  represented  by  a  class  of  integrals.  The  most  useful 
information  for  our  purpose  (understanding  the  behavior  of  the  difference  approximations) 
is  contained  in  graphs  of  these  integrals,  since  the  solution  of  a  difference  approximation 
is  closely  modeled  by  these  integrals. 

We  will  first  transform  the  general  canonical  form  given  in  section  4  to  both  bring  out 
the  independent  parameters  and  to  place  it  in  a  form  with  convenient  limits.  Specifically, 
we  would  like  to  show  the  scalar  wave  equation  limit  as  a  way  of  showing  the  similarities 
in  the  behavior  of  approximations  to  the  two  problems. 

The  canonical  form  integral  in  section  4  is 


yft/K 

J-*/ h 


(  ’exp 


+  *'£(<*> 


1)  4-  ia£{ih£)p  + 


We  will  take  p  even  and  q  odd,  and  s  =  1  or  0.  After  a  little  bit  of  algebraic  manipulation, 
we  may  write  this  as 


+  ixy  + 


iyp+ 1 
P  +  1 


ayq+ 1 
9  +  1 


Ay. 


(The  limits  of  integration  have  been  extended  to  ±oo.)  The  values  for  a,  0,  and  x  are 


a  = 


x  = 


\)wa  P+7 


62 


where 


a  =  Ip  +  i)a(~1)5 
b  =  -(</  +  l)t(-l)V. 

(u  corresponds  to  c  and  b  corresponds  to  t  in  Chin  and  Hedstrom  [1978].)  Note  that  ft 
is  proportional  to  the  square  of  the  cell  wave  coupling  number  H  discussed  in  section  4, 
and  that  only  ft  contains  any  contribution  from  p.  Thus  the  coupling  in  the  telegrapher's 
equation  shows  up  only  in  the  ft  term.  Further,  we  have  dropped  leading  constants.  We 
define  U -,(n,  ft,p,q-,  z)  as  this  integral.  The  advantage  of  this  form  is  that  I,  .,((»,  0,  p,  </-,  x) 
is  equal  to  the  generalized  Airy  functions  tabulated  by  Chit,  and  Hedstrom  [1978]  for  the 
scalar  wave  equation  (our  p  is  their  p—  1,  our  q  is  their  q  —  I).  This  will  allow  us  to  both 
check  our  graphs  against  their  published  solutions  and  to  see  how  the  ft  term  changes  the 
behavior  of  ft,  p,  q ;  z). 

The  graphs  are  organized  as  follows.  Each  page  has  graphs  for  fixed  a,  p  and  q. 
The  left  column  shows  l0(o,  ft,  p,  q-t  x)  for  ft  =  0,  0.01,  and  0.0f>.  The  right  column  gives 
li  («,  ft,  p,  q;  x)  for  the  same  values  of  ft-  The  left  column  represents  the  leading  term  in  the 
asymptotic  expansion  for  the  difference  solution. 

The  last  two  pages  of  these  graphs  show  the  special  case  where  there  is  no  ;>  term 
These  represent  the  odd  order  methods  (such  as  upstream  differencing). 

To  use  these  graphs,  the  following  should  be  kept  in  mind.  First,  in  comparing  the  lo 
graphs  with  the  computed  solutions,  the  left  column  should  be  multiplied  by  -l  (flipped 
over)  since  the  integral  for  the  solution  in  sections  3  and  4  is  the  negative  ot  the  form  used 
here.  Second,  except  for  the  last  two  pages  (the  ‘p-term  missing'  graphs),  these  graphs 
are  appropriate  for  even  order  methods  with  artificial  viscosity  (or  no  viscosity  if  a  —  0). 
Increasing  values  of  a  represent  increasing  artificial  viscosity. 
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Explanation  and  Use  of  the  Graphs 


These  graphs  show  the  qualitative  similarity  of  the  behavior  of  difference  approxima¬ 
tions  to  the  scalar  wave  equation  (represented  by  0  =  0)  and  the  telegrapher’s  equation 
(0  7^  0)  near  a  discontinuity. 

The  first  thing  to  notice  from  these  graphs  is  the  qualitative  similarity  for  graphs  that 
differ  only  in  0.  The  only  noticeable  change  in  the  graphs  as  0  increases  from  0  is,  for 
l0,  a  change  in  slope  on  the  left  of  the  jump  and,  for  I,,  a  change  in  the  zero  level  of  the 
oscillations  on  the  left  of  the  jump.  From  this  we  can  conclude  that  near  the  discontinuity, 
oscillations  in  the  approximate  solution  to  (1.1)  are  qualitatively  similiar  to  those  of  the  wave 
equation  (0  =  0). 

We  may  also  compare  these  graphs  to  the  calculations  of  actual  methods  presented 
above.  For  example,  the  graph  for  p  —  2,  q  ~  3,  a  =  1/2,  and  0  0.05  and  the 

Lax-Wendroff  solution  are  very  similiar.  The  graph  for  p  =  2,  n  =  0,  and  0  =  0.05  shows 
the  long  trail  of  oscillations  which  we  see  in  both  Leapfrog  and  the  second  order  centered 
difference  approximation.  The  graph  for  p  =  4,  a  =  0,  and  0  =  0.05  for  nondissipative 
fourth  order  methods  shows  both  the  long  trail  of  oscillations  and  the  small  dip  ahead  of 
the  front  which  we  observed  in  the  graph  of  the  solution  by  the  fourth  order  method. 

As  an  illustration  of  the  use  of  these  tables,  we  can  see  that  because  the  similiarity  in 
the  behavior  of  the  approximate  solutions  is  so  strong,  the  same  artificial  viscosities  should 
work  for  the  telegrapher's  equation  as  for  the  scalar  wave  equation  (in  terms  of  improving 
the  solution  near  the  discontinuity). 


Application  to  Mesh  Refinement 

We  present  here  an  application  of  the  results  in  section  4.  We  will  apply  mesh 
refinement  to  the  example  used  by  Majda  and  Osher  [1977],  Their  problem  was 


() 

Ot 


°)UUi)+(  o 


ui(z,0)  — 


if  z  <  0 
if  z  >  0 


u2(x,  0)  =  0. 


(5.5) 


We  will  use  a  fourth  order  method  recommended  by  Apelkrans  [1969]  in  the  region 
away  from  the  discontinuities,  and  the  second  order  Lax-Wendroff  method  around  the 
discontinuities  and  at  the  boundaries. 


The  form  of  mesh  refinement  that  we  use  is  described  in  detail  by  Berger,  Gropp, 
and  Oliger  [unpublished];  the  mesh  refinement  program  is  due  to  Berger  [unpublished]. 
Basically,  this  form  of  mesh  refinement  places  refined  grids  dynamically  over  the  problem 
domain.  The  grids  may  be  refined  to  an  arbitrary  degree. 

Since  we  must  compute  on  a  finite  domain,  we  must  introduce  a  boundary  and  a 
boundary  approximation.  From  the  results  in  Gustafsson  [1975],  we  know  that  the  boundary 
approximation  must  be  at  least  third  order  accurate  if  we  want  fourth  order  accuracy  for 
the  computation  as  a  whole.  There  are  two  common  ways  to  achieve  this:  use  a  third 
order  uncentered  scheme  such  as  that  recommended  in  Oliger  [1974]  for  the  scalar  wave 
equation,  or  use  a  second  order  method  and  mesh  refinement,  as  recommended  in  Oliger 
[1976],  We  have  chosen  the  second  approach  as  it  is  the  easiest  to  apply  in  this  problem. 

Our  algorithm  for  this  problem  is  the  following:  we  apply  a  fourth  order,  dissipative 
scheme  on  a  grid  with  mesh  size  hc  except  near  the  boundaries  and  the  discontinuities 
at  x  =  ±t,  and  we  apply  a  second  order,  dissipative  method  on  a  grid  with  mesh  size 
hf  —  h2c  at  those  places. 

The  first  question  which  may  come  to  mind  is  the  computational  complexity  of  such 
an  algorithm.  We  will  first  estimate  the  total  number  of  mesh  points,  both  coarse  and  fine, 
needed  by  the  algorithm.  From  that  result,  we  will  be  able  to  compute  the  computational 
effort  required  by  this  algorithm.  There  are  three  contributions  to  the  number  of  mesh 
points:  the  coarse  grid,  the  fine  grid  at  the  boundaries,  and  the  fine  grid  around  the 
discontinuities.  The  coarse  grid  has  0(l//ic)  mesh  points.  The  fine  grid  at  the  be  ^daries 
has  0(\/hf)  x  \hc  —  0(1 /hc)  mesh  points,  since  the  boundary  grid  only  extends  over  4 
coarse  grid  points. 

To  determine  the  number  of  fine  grid  mesh  points  around  the  discontinuity,  we  must 
determine  the  width  of  the  front.  We  make  the  two  assumptions:  the  region  of  refinement 
is  asymptotically  of  the  same  size  as  the  width  of  the  front,  and  the  derivative  discontinuity 
at  x  =  t  is  asymptotically  no  larger  than  the  discontinuity  at  r.  —  —t.  On  the  fine  grid,  the 
results  of  section  4  show  that  the  width  of  the  front  is  0(1 i^/3).  But  since  hf  —  h~,  the 
width  of  the  front  on  the  scale  of  the  coarse  grid  is  0(hy*).  Thus,  as  /.  0,  the  number 

of  coarse  mesh  points  that  the  front  is  spread  across  tends  to  zero.  Because  of  the  way 
that  the  mesh  refinement  program  works,  there  is  a  minumum  of  seven  coarse  grid  points 
(six  cells)  which  must  be  refined  around  the  fronts.  Thus,  the  number  of  fine  grids  points 

87 


around  the  discontinuities  is  0 ( I / /t / )  X  f>hc  —  0(l//»c).  Combining  these  gives  0(l//ir) 
total  mesh  points. 

The  computational  work  is  easily  derived  from  the  above  analysis.  To  integrate  to  time 
To  with  time  step  kc  on  the  coarse  grid  and  time  step  kf  ~  k'f  on  the  fine  grid,  we  need 


0(1/*.)  Xji 

0(1/*.)  xg 


on  the  coarse  grid 
on  the  fine  grid 


or  0(\/h 2)  +  0(1/ h*)  work  (assuming  kc/hc  —  constant).  Thus  we  see  that  the  work  on 
the  fine  grid  dominates  the  work  on  the  coarse  grid.  These  estimates  are  observed  in 
computations. 

These  results  show  that  this  algorithm  is  practical  in  space  (it  uses  no  more  space 
than  the  coarse  grid  would  alone),  but  less  practical  in  time.  The  advantage  of  this  method 
is  the  ability  to  compute  a  more  accurate  solution  for  a  given  coarse  grid  than  is  possible 
without  refinement. 

For  our  computational  example,  we  computed  the  solution  to  (5.5)  at  t  =  0.5  with 
hc  —  0.02  and  hc  =  0.01.  The  time  step  fc  was  taken  as  hc/ 2.  Computations  with  and 
without  mesh  refinement  at  the  discontinuities  were  done;  their  results  are  presented  in  the 
following  table.  These  choices  of  t  and  h  were  made  to  facilitate  comparisons  with  the 
results  in  Majda  and  Osher  [1977], 
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X 

solution 

U\ 

MR: 

h  -  0.02 

MR: 

h  =  0.01 

no  MR: 
h  =  0.02 

no  MR: 
h  =  0.01 

~0.8  ~ 

8.775826E-1 

3.9E-7 

<  IE-7 

3.0E-7 

<1E-7 

07 

8.775826E-1 

C1E-7 

<  IE-7 

-3.1E-5 

<  IE-7 

0.6 

8.775826E-1 

2.5E-6 

<  IE-7 

-3.8E-3 

2. IE-5 

-0.5 

-1.224174E-1 

-6.6E-1 

-6.6E-1 

-5.7E-1 

5.7E-1 

_  _ 

0.4 

-9.879886E-2 

-3.8E-3 

3.5E-5 

-4.2E-2 

2.7E-2 

~  -0.3_ 

-  7.78861 3E-2^ 

7.9E-5 

-7. IE-6 

T1E2 

1 .2E-3 

T).2 

-5.956894E-2 

8.0E-5 

-1.9E-7 

-2.4E-3 

-6.2E-5 

-01 

-4.377426E-2 

-6.6E-7 

-2.6E-7 

1.9E-4 

aoE  6  ’ 

0.0 

-3.044362E-2 

2.8E-6 

3.2E-7 

9.9E  6 

2.0E-6 

0.1 

-1.953720E-2 

-4.9E  7 

4.9E-7 

7.6E-6 

2.0E-6 

0.2 

-1.103355E-2 

-1.5E-6 

5.1E-7 

8.1E-6 

2.0E-6 

0.3 

-4.929534E-3 

-6.0E-7 

2.6E-8 

8.4E-6 

2.0E-6 

0.4 

-1.234040E-3 

-3.1E-9 

2.7E-8 

8.6E-6 

1 ,5E-6^ 

0.5 

0.0 

_ 

4. IE-7 

6.6E-8 

1  3E  5 

3.9E-6 

0.6 

0.0 

1.0E-11 

-5.2E-16 

1.5E-7 

-1.5E-9 

0.7 

6.0 

-  -  .  -  .  . 

8.0E-15 

-3.9E-21 

4.1  E- 1 1 

6.0E  1 5 

0.8 

0.0 

4.3E-18 

-3.3E-28 

3.2E-13 

-1 .9E-20 

Table  4:  Errors  in  the  computed  solution  at  x  =  -.8(1). 8  for 
the  fourth  order  method,  with  and  without  mesh  refinement.  The 
columns  labeled  MR  are  the  results  of  the  mesh  refinement  com¬ 
putations. 

It  is  clear  from  this  table  that  mesh  refinement  significantly  improves  the  accuracy  of 
the  results.  Further,  the  rate  of  convergence  is  higher  for  the  mesh  refinement  calculation 
than  for  the  unrefined  calculation. 

This  method  is  not  the  most  efficient  way  to  solve  (5.5).  (For  this  problem,  front 
tracking  is  more  efficient.)  It  does,  however,  illustrate  the  use  of  mesh  refinement,  coupled 
with  knowledge  of  the  behavior  of  the  problem,  to  get  a  high  accuracy  solution. 
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6.  Conclusion 


Summary  of  Results 

We  have  shown  that  the  asymptotic  width  of  the  front  is  hp/l(pi',\  where  p  is  the 
order  of  the  approximation,  for  a  fairly  general  class  of  difference  approximations  to  the 
telegrapher’s  equation.  This  behavior  is  very  similiar  to  the  behavior  near  the  front  for  the 
scalar  wave  equation,  and  shows  that  the  effects  of  numerical  dispersion  and/or  dissipation 
near  the  discontinuity  overwhelm  any  other  numerical  artifacts  which  may  be  present,  such 
as  the  effect  discovered  by  Mnida  and  Osher  [1977], 

The  canonical  forms  that  we  constructed  contain  a  great  deal  of  information  on  the 
detailed  behavior  of  the  solution.  Given  a  difference  scheme  that  fits  the  framework  in 
section  4,  we  can  use  the  graphs  of  the  canonical  forms  dislayed  in  section  5  to  get  an 
idea  of  how  the  solution  to  the  difference  scheme  will  behave  near  a  discontinuity. 

The  detailed  examination  of  the  non  dissipative  second  order  difference  scheme  studied 
in  section  3  shows  that  the  canonical  forms  displayed  in  section  5  are  actually  a  good 
approximation  to  the  exact  canonical  form  near  the  front  (see  figures  10  and  11).  These 
results  also  show  that  the  model  equation  is  a  good  approximation  near  the  front  (because 
h£  for  the  saddle  points  is  <>(  I )  near  the  front).  Finally,  the  canonical  form  for  the 
approximation  in  section  3  shows  that  away  from  the  front,  the  model  equation  formulation 
is  not  accurate  (figure  11). 

Another  difficulty  in  determining  the  solution  away  from  the  front  is  the  complicated 
behavior  of  the  unexponentia'ed  term  in  the  integral  for  the  solution.  Even  for  the  exact 
(continuous)  problem,  this  is  difficult  to  deal  with. 

In  section  4  we  also  investigated  the  effect  of  the  coupling  term.  We  showed  that 
the  effect  of  the  coupling  term  is  controlled  by  the  cell  wave  coupling  number,  which  is 
equal  to  the  product  of  the  step  size  and  the  coupling  strength,  divided  by  the  speed  of 
propagation  (hp/r).  As  long  as  il  is  small,  the  behavior  near  the  front  of  approximations  to 
the  telegrapher's  equation  is  similar  to  the  behavior  of  approximations  to  the  scalar  wave 
equation.  This  result  is  remarkable  because  the  results  of  Majda  and  Osher  [1977)  show 
that  away  from  the  front,  the  behavior  of  approximations  to  these  two  equations  is  very 
different.  The  analysis  in  section  4  showed  further  that  the  effect  of  the  coupling  term  is 
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noticeable  only  very  near  the  front  (where  it  affects  the  I1  or  forward  moving  term)  and 
well  away  from  the  front  (where  the  noise  generated  at  the  front,  moving  away  from  the 
front,  shows  up  in  the  /  term). 

An  interesting  feature  of  difference  approximations  to  the  telegrapher  s  problem  which 
we  have  shown  is  the  different  frequency  content  of  the  solution  of  the  differential  equation 
as  compared  to  the  solution  of  difference  approximations.  For  the  differential  equation 
saw  in  section  2  that  the  leading  terms  in  the  expansion  of  the  inverse  Fourier  transform 
were  at  infinite  frequency.  Specifically,  the  saddle  points  were  converging  on  infinity  m  the 
frequency  (£)  plane  as  the  front  was  approached.  In  contrast,  we  saw  in  sections  3  and  4 
that  the  saddle  points  for  the  difference  approximations  converged  on  zero  in  the  vp  piano 
It  is  this  fact  that  makes  the  mode!  equations  analysis  in  section  4  work,  since  the  mo  !>4 
equation  is  essentially  an  expansion  in 

The  computations  presented  in  section  5  support  our  theory.  In  every  case  the 
asymptotic  width  of  the  front  matched  the  hv/{p 1  prediction.  The  graphs  of  the  canonic  al 
forms  for  various  difference  schemes  matched  the  qualitative  behavior  near  the  front  for 
those  difference  schemes. 

Implications 

The  results  in  this  thesis  show  that  the  effect  of  numerical  dispersion  and  hs*.  pat  a 
near  a  discontinuity  overwhelms  the  numerical  artifact  introduced  by  coupling  r  " 
equations.  This  implies  that  the  effect  discovered  by  Majda  and  Oshei  ina,  be  ivgk-- 
if  (a)  accurate  computation  of  the  front  is  of  primary  importance  and  (b)  second  )-_»» 
accuracy  is  sufficient  in  the  smooth  regions.  In  this  case,  mesh  refinement  y..u!od  L,  vne 
knowledge  of  the  width  of  the  front,  can  be  used  to  accurately  compute  the  discGr.teuut , 
Further,  artificial  viscosities  which  are  good  for  the  scalar  wave  equation  will  be  good  hne 
since  artificial  viscosities  are  designed  to  reduce  certain  high  frequencies  and  m.r  n  .  . d ! ■ 
show  that  the  frequency  of  the  oscillations  near  the  front  are  nearly  the  same  f,v  the  m  ,0  o 
wave  equation  and  the  telegrapher's  equation. 

If  an  accurate  solution  is  needed  everywhere,  mesh  refinement  can  bo  c  ed  as 
above,  with  the  added  requirement  that  the  solution  near  any  discontinuity  bo  calculated 
accurately  enough  to  prevent  the  contamination  cf  the  smooth  pait  of  the  solut.on  I  ‘a-- 
war,  demonstrated  in  section  5. 
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Topics  for  Further  Study 

The  study  of  even  such  a  narrow  topic  as  this  is  never  complete;  below  are  a  number 
of  questions  that  arose  during  the  research  into  this  problem. 

Multiple  time  scales  —  Many  transport  equations  that  arise  in  practice  involve  two 
radically  different  materials.  For  example,  astrophysical  transport  problems  often  involve 
the  interaction  of  electromagnetic  radiation  with  gas.  The  scale  of  time  during  which  the 
radiation  moves  is  orders  of  magnitute  different  from  the  scale  of  time  on  which  the  gas 
moves.  A  current  approach  to  these  problems  is  to  solve  seperate  equations  for  the  gas 
and  the  radiation,  then  coupling  them  together  at  long  (for  the  radiation  intervals).  Another 
approach  to  these  problems  is  to  use  stiff-integrators  to  integrate  the  equations  in  time. 
Questions  that  arise  here  are  (a)  how  accurate  is  this?  (b)  how  do  these  algorithms  behave 
in  the  presence  of  discontinuities?  (c)  are  there  better  ways? 

Asymptotic  solution  away  from  the  front  —  For  the  general  problem  we  were  unable  to 
find  an  asymptotic  representation  for  the  solution  except  near  the  discontinuity.  The  model 
equation  approach  does  not  seem  powerful  enough  to  ever  provide  a  way  to  represent 
the  solution  away  from  the  front.  Is  there  a  general  technique?  (Specific  methods  can  be 
studied,  as  we  did  in  section  3,  and  some  specific  problems  can  be  studied  as  in  Majda 
and  Osher  [1977].) 

Initial  Boundary  Value  Problem  —  What  is  the  effect  of  coupling  through  boundary 
conditions?  The  methods  used  in  this  thesis,  with  Fourier  transforms  replaced  by  Laplace 
transforms,  would  give  an  answer. 
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Appendix  A:  Two  Integral  Identities 


We  prove  in  this  appendix  two  integral  identities  that  we  used  in  section  2.  We  first 
show  that 


0(0  exp  =  0 


(A.1) 


where 


and 


0(0  = 


1  +  ty/4f8  -  1  -  2t£ 
~  1 


^0  =  ^V^rrT  +  *«0 


As  always,  the  integral  is  the  Cauchy  principle  value,  and  we  pass  above  the  branch  cut 
[—1/2, 1/2].  To  make  the  proof  somewhat  clearer,  we  will  assume  that  t  =  1.  It  will  be 
clear  that  the  proof  will  hold  for  any  t  >  0. 

The  proof  is  a  simple  exercise  in  using  the  Cauchy  integral  theorem.  We  close  the 
contour  in  the  upper  half  plane  with  a  semicircle  of  radius  R.  The  proof  consists  of  showing 
that  as  R  -*  oo,  the  contribution  to  the  integral  over  the  closed  contour  from  the  semicircle 
tends  to  zero.  Finally,  it  is  demonstrated  that,  despite  appearances,  g(£)  is  regular  at  the 
origin  and  only  of  the  order  of  (£  ±  |)_I/2  at  £  and 

Along  the  semicircle,  we  have  £  =  Re'9.  We  will  make  this  substitution  in  the  integral 
as  we  go  along  in  the  proof. 

We  start  by  finding  a  bound  on  K(  j>/4£2  —  1).  This  gives  us  the  magnitude  of  the 
exponential  in  the  integral.  First,  we  see  that 


where  3R(z)  is  the  real  part  of  z  and  3(z)  is  the  imaginary  part  of  z.  This  is  equal  to 


— R  sin  0  8? 


(\A  -  1ST  )  - ' R ' cm0  3(\/1  -  ^  )• 


Now,  a  quick  look  at  Ahlfors  [1979]  shows  us  that 
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*  \/l- 


I  —  E2t3i  .  ft  4.  _J _ c°‘  21 

1  4R1  +  V  T  >««4  a*1 


9  V  1  — 


£SJ_21  _  I  4.  .  /I  4.  — 1 _ col  2« 

4«l  1  r  V  1  ^  i«ft4  a*’ 


We  bound  1  -  J  from  below.  The  argument  of  the  inner  square  root  satisfies 

(  1  V  1  cos  29  /  I  y 

V  +  4«V  -  +  16ft4  ~  2ft2  -  (  _  4K2] 


so  that 


«  2‘*  ^  -  ,  /*  C4«>*  +  1  4fla 


"(V'-w  )* 


for  ft  >  1.  Similiarly,  we  get  a  bound  on  the  imaginary  part 


for  ft  >  5.  Thus, 


<p(0U))l  <  exp^-ftsinfl  +  ~ exp(t'wC)l  <  2exp(~f  +  w)sintfj.  (A.3) 


We  also  need  a  bound  on  |g(£)|. 


|g(fte*)|  < 

We  can  bound  Iv^2  -  t|  =  y/\iR*&*  -  1|  =  |16ft4  -  8fta  cos  20  +  1|*  by 


ZR  >  J\\R*e2it  -  1|  >  ft 


for  ft  >  t.  Thus, 


for  ft  >  1. 


k(Rc‘#)l  <  | 


Now,  using  (A.3)  and  (A.4),  we  can  bound  the  contribution  to  the  integral  over  the 
semicircle  by 


8  f  e-,*(i+w)!,inJ  dO 
J  0 
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We  can  bound  this  integral  by  dividing  it  into  two  parts.  The  first  part  is  0  <  0  <  c/16. 
The  part  of  the  integral  over  this  range  erf  9  is  no  greater  than  c/2  for  R  >  1.  The  second 
part  is  for  c/16  <  0  <  */2.  In  this  range  of  0,  we  can  pick  R  to  satisfy 

Then  the  integrand  is  less  than  c/8w,  and  the  contribution  to  the  integral  is  less  than  c/2. 
Thus,  by  choosing  R  large  enough,  we  can  make  the  integral  (A.5)  as  small  as  we  want 
(for  u)  >  -1/2). 

Next,  we  must  treat  the  singularities  of  9(C).  By  inspection,  it  is  clear  that  there  are 
at  most  three  points  of  singularity  in  g(().  We  can  investigate  these  by  using  macsyma 
(Mathlab  Group  [1975])  to  expand  g{()  about  the  points  1,  J,  and  We  have 

9(0*-j~y-C*  +  —  for  £  »  0 

,_1  ,•  (5-V?±T 

g(£)=±  — - = - +  •••  for  (  «  T j 

^ T  ‘ 

Thus  there  is  no  contribution  to  the  Cauchy  principle  value  of  (A.1)  from  any  of  these,  as 
they  are  weaker  than  simple  poles.  The  relation  (A.1)  is  now  proved. 

A  stronger  result  can  be  proved  from  (A.3)  and  (A.2)  by  noticing  that  the  5  in  ( £  +  w) 
is  really  the  (1  -  1  /(2/?a))l/a  in  (A.2).  Therefore,  as  R  -+  00,  the  bound  on  (A.3)  is 
asymptotically 

|cxp(*(0)i  &  2cxp(-K(l  +u)s\nO) 

for  w  >  - 1 .  This  matches  what  we  would  expect  to  get  from  an  asymptotic  expansion  of 
<f>,  since  <j>(()  ~  *(1  +  w)£  for  |£|  large. 


We  next  prove  that 


dx 


=  2xt(-l)^Q  '  1,(2 sffiS 


(A.6) 


for  integer  v  >  0  and  and  p  and  7  pure  imaginary.  Of  course,  we  show  this  in  the 
generalized  sense  of  Gel' land  and  Shilov  [1964]. 

We  start  by  assuming 


*-,e-1 r*  =  2iri\i)(2y/p^). 


(A.7) 
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We  proceed  to  show  that  we  can  get  (A.6)  from  this  by  formally  differentiating  each  side  of 
this  with  respect  to  7,  passing  the  derivative  through  the  integral. 

First,  we  note  that 


Xu-ie-jx-fi/x  d  =(_!)*  I  x-i  ^Le-r *-P/x  dx 

1  '  y_oo  dY 

JV  f  +  °° 

v  1  dY  J- 00 

=(-l),'2irt^:lo(2V^r). 

We  can  evaluate  the  derivative  by  an  induction  argument. 

From  Abramowitz  and  Stegun  [1972],  we  have  the  formula 


{zdz){x^[z))~ 


or 


v  ,  .  ,  1  d,  .  .  1,  .. 

+  ZZ+TTZ^^*)  ~~  ZZ+T^+d*)- 


gW*  -  •  '  zv+\  dz  V  ’  zu+i 
Now,  let  z  =  2'v/^7,  and  multiply  the  above  by  (2/J)*'+1.  This  leaves  us  with 

+ (;)+ = <«> 


2\Z?T 

Now,  we  start  with  the  induction  proof.  The  first  term,  u  =  1,  we  .can  get  from  the 
fact  that  Iq  =  Ij.  This  gives  us 

which  is  the  first  term  in  the  induction.  For  the  nth  term  in  the  induction,  we  need  to  show 

i(02flU-1(2^)  =  (0’w2v^). 

We  do  this  by  performing  the  indicated  differentiation  and  then  applying  (A.8). 

=(^)’u(2vW). 


Here  we  have  substituted  n  -  1  for  v  in  (A.8). 
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We  need  only  prove  (A.7)  to  prove  (A.6).  We  make  the  change  of  variable  x  =  iu 
(A.7).  This  gives  us 


in 


From  page  245  of  Erdelyi  [1954],  this  is  just  an  inverse  Laplace  transform  with  value 
2iri\Q(2\/jFi).  Note  that  for  p  and  =  0,  a  singular  functional  which  coincides  with  zero 
for  p  and  7  0  must  be  added  to  (A.6).  This  completes  our  proof  of  (A.6). 
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Appendix  B:  Computing  the  Exact  Solution 


This  appendix  discusses  the  method  used  to  compute  the  exact  solution  to  the 
continous  model  problem  in  section  2.  From  section  2,  we  know  that  we  need  to  evaluate 
the  integral 

1  r+°° 

u_  =-^-  /  a_(l  +  d£ 

J  —  oo 

2ir  J_oo  4$s/4^2  - 1  4 

We  need  to  perform  this  integration  numerically.  First,  we  will  deform  the  path  of  integration 
to  a  path  “close”  to  the  path  of  steepest  descent.  The  figures  below  show  the  qualitative 
behavior  of  the  path  of  steepest  descent;  we  will  use  these  to  guide  us  in  our  choice  of 


From  these  figures  it  is  dear  that  the  path  to  take  is  a  circle  about  the  origin.  In 
deforming  the  path  to  this  circle,  we  will  pick  up  —  {  from  the  pole  at  the  origin.  We 
can  show  by  techniques  similiar  to  those  in  Appendix  A  that  there  is  no  contribution  from 
closing  the  contour  in  the  lower  half  plane. 

The  resulting  integral  is 


±  r  i 

2rr  Jo 


2*  Jo  4>/4£2  -7 


dO 


(8.1) 


where  £  =  relt  and  r  has  been  chosen  as  max(w/2\A  —  u*,  |  +  e),  c  some  small  positive 
number.  This  choice  of  r  makes  the  path  pass  through  the  saddle  points  unless  the  path 
would  pass  within  e  of  the  branch  cut  [-£,  J],  in  which  case  the  path  is  expanded. 

This  integral  has  a  special  form.  The  integrand  is  periodic  analytic  in  a  strip  around 
the  region  of  integration.  It  is  well  known  that  the  trapezoidal  rule  converges  very  quickly 
for  periodic  functions  (Davis  and  Rabinowitz  [1975]);  in  particular,  for  periodic  analytic 


functions  which  satisfy  certain  technical  assumptions,  the  convergence  is  exponential  (see 
Rabinowitz  [1968]  and  Lyness  and  Delves  [1967]).  Using  this  fact,  the  integral  (B.1)  was 
integrated  with  the  trapezoidal  rule. 

We  would  like  to  point  out  a  few  features.  Along  the  steepest  descent  path,  the 
imaginary  part  of  <£_  is  constant.  This  means  that  the  integrand  should  not  have  serious 
oscillations  along  the  path  of  integration.  Secondly,  since  most  of  the  effect  is  concentrated 
at  the  saddle  points,  we  would  expect  an  automatic  quadrature  routine  which  places  points 
adaptively  not  to  be  that  much  less  efficient  than  the  trapezoidal  rule.  Experimentally  this 
is  the  case.  The  routine  QAQS  (de  Doncker  [1978])  was  roughly  as  fast  as  a  simple 
trapezoidal  sum.  Finally,  the  choice  of  e,  the  minimum  distance  from  the  singularities,  can 
have  a  big  effect.  Since  there  are  singularities  at  £  =  ±5,  we  should  keep  the  path  well 
away  from  these  points.  A  choice  of  j  for  t  seems  to  work  well. 
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